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 : 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 :
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 46 : SUBROUTINE create_A(fm_S_ia, fm_S_bar_ij, fm_S_ab, &
97 46 : fm_A, Eigenval, unit_nr, &
98 46 : 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 46 : INTEGER, ALLOCATABLE, DIMENSION(:) :: eig_offsets, n_ov, offsets
113 : INTEGER, DIMENSION(4) :: reordering
114 46 : 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 46 : CALL timeset(routineN, handle)
125 :
126 46 : nspins = SIZE(homo)
127 230 : ALLOCATE (n_ov(nspins), offsets(nspins), eig_offsets(nspins))
128 46 : 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 46 : eig_offsets(1) = 0
131 58 : DO isp = 2, nspins
132 58 : eig_offsets(isp) = eig_offsets(isp - 1) + homo(isp - 1) + virtual(isp - 1)
133 : END DO
134 :
135 46 : NULLIFY (dft_control, tddfpt_control)
136 46 : dft_control => bse_env%dft_control
137 46 : tddfpt_control => dft_control%tddfpt2_control
138 :
139 46 : 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 92 : SELECT CASE (bse_env%bse_spin_config)
145 : CASE (bse_singlet)
146 46 : alpha = 2.0_dp
147 : CASE (bse_triplet)
148 46 : alpha = 0.0_dp
149 : END SELECT
150 : ! For open-shell (nspins>1): each spin block contributes once; SPIN_CONFIG is ignored.
151 46 : IF (nspins > 1) THEN
152 : CALL cp_warn(__LOCATION__, &
153 12 : "BSE: SPIN_CONFIG ignored for open-shell reference; using alpha=1.")
154 12 : alpha = 1.0_dp
155 : END IF
156 :
157 46 : IF (bse_env%screening_method == bse_screening_alpha) THEN
158 2 : alpha_screening = bse_env%screening_factor
159 : ELSE
160 44 : 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 46 : NULLIFY (blacs_env)
165 46 : 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 46 : para_env=fm_S_ia(1)%matrix_struct%para_env)
175 46 : CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb")
176 46 : 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 46 : 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 46 : IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN
186 42 : 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 12 : para_env=fm_S_ia(1)%matrix_struct%para_env)
192 12 : CALL cp_fm_create(fm_S_joint, fm_struct_S_joint, name="fm_S_ia_joint")
193 12 : CALL cp_fm_set_all(fm_S_joint, 0.0_dp)
194 12 : 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 12 : matrix_c=fm_A)
198 12 : CALL cp_fm_release(fm_S_joint)
199 12 : 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 30 : beta=0.0_dp, matrix_c=fm_A)
205 : END IF
206 : END IF
207 :
208 46 : 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 104 : DO isp = 1, nspins
216 104 : 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 56 : para_env=fm_S_ab(isp)%matrix_struct%para_env)
220 56 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab")
221 56 : 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 56 : beta=0.0_dp, matrix_c=fm_W)
227 56 : 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 56 : row_offset=offsets(isp), col_offset=offsets(isp))
231 56 : 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 2 : 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 32 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
238 : tddfpt_control%do_bse_gw_only) THEN
239 4 : NULLIFY (ex_env)
240 4 : ex_env => bse_env%exstate_env
241 4 : IF (.NOT. tddfpt_control%do_bse_gw_only) THEN
242 12 : ALLOCATE (ex_env%bse_w_matrix_MO(1, 1))
243 12 : ALLOCATE (ex_env%bse_a_matrix_MO(1, 1))
244 4 : CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_W)
245 4 : CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_A)
246 4 : CALL cp_fm_to_fm(fm_W, ex_env%bse_w_matrix_MO(1, 1))
247 4 : IF (tddfpt_control%do_bse_w_only) THEN
248 2 : CALL cp_fm_to_fm(fm_A_copy, ex_env%bse_a_matrix_MO(1, 1))
249 : ELSE
250 2 : 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 56 : CALL cp_fm_release(fm_W)
256 112 : CALL cp_fm_struct_release(fm_struct_W)
257 : END IF
258 : END DO
259 :
260 46 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
261 3 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
262 : END IF
263 46 : 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 46 : 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 46 : IF (.NOT. tddfpt_control%do_bse) THEN
271 6440 : DO ii = 1, nrow_local_A
272 6396 : i_row_global = row_indices_A(ii)
273 2110312 : DO jj = 1, ncol_local_A
274 2103872 : j_col_global = col_indices_A(jj)
275 2110268 : 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 4826 : isp = nspins
278 4826 : DO k_isp = 1, nspins - 1
279 4826 : 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 3596 : k_ov = i_row_global - offsets(isp)
285 3596 : i_occ_row = occ_of_ia(k_ov, virtual(isp))
286 3596 : a_virt_row = virt_of_ia(k_ov, virtual(isp))
287 : eigen_diff = Eigenval(eig_offsets(isp) + a_virt_row + homo(isp)) - &
288 3596 : Eigenval(eig_offsets(isp) + i_occ_row)
289 3596 : 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 46 : IF (nspins == 1) THEN
297 34 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
298 : tddfpt_control%do_bse_gw_only) THEN
299 4 : sizeeigen = SIZE(Eigenval)
300 12 : ALLOCATE (ex_env%gw_eigen(sizeeigen))
301 126 : ex_env%gw_eigen(:) = Eigenval(:)
302 : END IF
303 : END IF
304 :
305 46 : CALL cp_fm_struct_release(fm_struct_A)
306 46 : DEALLOCATE (n_ov, offsets, eig_offsets)
307 :
308 46 : CALL cp_blacs_env_release(blacs_env)
309 :
310 46 : CALL timestop(handle)
311 :
312 138 : 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 22 : SUBROUTINE create_B(fm_S_ia, fm_S_bar_ia, fm_B, &
330 22 : 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 22 : CALL timeset(routineN, handle)
349 :
350 22 : nspins = SIZE(homo)
351 88 : ALLOCATE (n_ov(nspins), offsets(nspins))
352 22 : CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
353 :
354 22 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
355 2 : 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 44 : SELECT CASE (bse_env%bse_spin_config)
361 : CASE (bse_singlet)
362 22 : alpha = 2.0_dp
363 : CASE (bse_triplet)
364 22 : alpha = 0.0_dp
365 : END SELECT
366 22 : IF (nspins > 1) alpha = 1.0_dp
367 :
368 22 : IF (bse_env%screening_method == bse_screening_alpha) THEN
369 2 : alpha_screening = bse_env%screening_factor
370 : ELSE
371 20 : 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 22 : 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 22 : para_env=fm_S_ia(1)%matrix_struct%para_env)
379 22 : CALL cp_fm_create(fm_B, fm_struct_B, name="fm_B_iajb")
380 22 : 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 22 : IF (nspins > 1) THEN
384 4 : 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 4 : para_env=fm_S_ia(1)%matrix_struct%para_env)
389 4 : CALL cp_fm_create(fm_S_joint, fm_struct_S_joint, name="fm_S_ia_joint")
390 4 : CALL cp_fm_set_all(fm_S_joint, 0.0_dp)
391 4 : 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 4 : matrix_c=fm_B)
395 4 : CALL cp_fm_release(fm_S_joint)
396 4 : 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 18 : beta=0.0_dp, matrix_c=fm_B)
402 : END IF
403 :
404 22 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
405 2 : 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 22 : IF (bse_env%screening_method /= bse_screening_rpa) THEN
411 44 : DO isp = 1, nspins
412 24 : 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 24 : para_env=fm_S_ia(isp)%matrix_struct%para_env)
418 24 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ibaj")
419 24 : 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 24 : beta=0.0_dp, matrix_c=fm_W)
424 24 : 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 24 : row_offset=offsets(isp), col_offset=offsets(isp))
428 24 : CALL cp_fm_release(fm_W)
429 68 : CALL cp_fm_struct_release(fm_struct_W)
430 : END DO
431 : END IF
432 :
433 22 : CALL cp_fm_struct_release(fm_struct_B)
434 22 : DEALLOCATE (n_ov, offsets)
435 :
436 22 : CALL timestop(handle)
437 :
438 44 : 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 46 : SUBROUTINE create_A_and_B(fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, &
459 92 : fm_S_bar_ij, fm_A, fm_B, do_abba, Eigenval, unit_nr, &
460 46 : 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 46 : 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 4 : homo, virtual, dimen_RI, bse_env)
476 4 : IF (do_abba) THEN
477 : CALL create_B(fm_S_ia, fm_S_ia, fm_B, homo, virtual, dimen_RI, unit_nr, &
478 4 : 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 42 : homo, virtual, dimen_RI, bse_env)
483 42 : IF (do_abba) THEN
484 : CALL create_B(fm_S_ia, fm_S_bar_ia, fm_B, homo, virtual, dimen_RI, unit_nr, &
485 18 : bse_env)
486 : END IF
487 : END IF
488 :
489 46 : 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 198 : 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 22 : CALL timeset(routineN, handle)
524 :
525 22 : IF (unit_nr > 0) THEN
526 11 : WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
527 22 : ' 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 22 : CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct)
541 22 : CALL cp_fm_to_fm(fm_A, fm_A_plus_B)
542 22 : CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct)
543 22 : CALL cp_fm_to_fm(fm_A, fm_A_minus_B)
544 22 : CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct)
545 22 : CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp)
546 22 : CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct)
547 22 : CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp)
548 :
549 22 : CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
550 :
551 22 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
552 2 : 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 22 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B)
557 22 : 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 22 : 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 22 : CALL cp_fm_create(fm_dummy, fm_A%matrix_struct)
567 : ! Create (A-B)^-0.5 (cf. EQ.II)
568 22 : 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 22 : CALL cp_fm_release(fm_dummy)
570 22 : IF (unit_nr > 0) THEN
571 11 : WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
572 22 : '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 22 : IF (eigvals_AB_diff(1) < 0) THEN
577 : CALL cp_abort(__LOCATION__, &
578 : "Matrix (A-B) is not positive definite. "// &
579 0 : "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 22 : 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 22 : 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 22 : fm_work_product)
592 :
593 : ! Release to save memory
594 22 : CALL cp_fm_release(fm_A_plus_B)
595 22 : CALL cp_fm_release(fm_A_minus_B)
596 :
597 : ! Now create full
598 22 : CALL cp_fm_create(fm_C, fm_A%matrix_struct)
599 22 : 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 22 : fm_C)
603 22 : CALL cp_fm_release(fm_work_product)
604 :
605 22 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
606 2 : 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 22 : CALL timestop(handle)
610 22 : 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 22 : SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
628 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
629 22 : 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 22 : 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 22 : CALL timeset(routineN, handle)
648 :
649 22 : nspins = SIZE(homo)
650 48 : n_ov_joint = SUM(homo*virtual)
651 :
652 22 : IF (unit_nr > 0) THEN
653 11 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
654 22 : '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 22 : CALL cp_fm_create(fm_eigvec_Z, fm_C%matrix_struct)
660 :
661 66 : ALLOCATE (Exc_ens(n_ov_joint))
662 :
663 22 : CALL choose_eigv_solver(fm_C, fm_eigvec_Z, Exc_ens, diag_info)
664 :
665 22 : IF (diag_info /= 0) THEN
666 : CALL cp_abort(__LOCATION__, &
667 0 : "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 22 : IF (Exc_ens(1) < 0) THEN
674 0 : 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 0 : "(A+B) is not positive definite.")
678 : END IF
679 : END IF
680 2650 : 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 22 : CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_C%matrix_struct)
692 22 : 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 22 : matrix_c=fm_mat_eigvec_transform_sum)
696 22 : CALL cp_fm_release(fm_sqrt_A_minus_B)
697 : ! This normalizes the eigenvectors
698 22 : 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 22 : CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_C%matrix_struct)
702 22 : 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 22 : matrix_c=fm_mat_eigvec_transform_diff)
706 22 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
707 22 : CALL cp_fm_release(fm_eigvec_Z)
708 :
709 : ! This normalizes the eigenvectors
710 22 : 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 22 : CALL cp_fm_create(fm_eigvec_X, fm_C%matrix_struct)
715 22 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_X)
716 22 : 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 22 : CALL cp_fm_create(fm_eigvec_Y, fm_C%matrix_struct)
721 22 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_Y)
722 22 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_Y, -1.0_dp, fm_mat_eigvec_transform_diff)
723 :
724 : !Cleanup
725 22 : CALL cp_fm_release(fm_mat_eigvec_transform_diff)
726 22 : CALL cp_fm_release(fm_mat_eigvec_transform_sum)
727 :
728 22 : 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 18 : .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 4 : .FALSE., mo_coeff, bse_env, unit_nr, fm_eigvec_Y)
736 : END IF
737 :
738 22 : DEALLOCATE (Exc_ens)
739 22 : CALL cp_fm_release(fm_eigvec_X)
740 22 : CALL cp_fm_release(fm_eigvec_Y)
741 :
742 22 : CALL timestop(handle)
743 :
744 154 : 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 22 : SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
758 22 : 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 22 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
771 : TYPE(cp_fm_type) :: fm_eigvec
772 :
773 22 : CALL timeset(routineN, handle)
774 :
775 22 : nspins = SIZE(homo)
776 52 : n_ov_joint = SUM(homo*virtual)
777 :
778 22 : IF (unit_nr > 0) THEN
779 11 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
780 22 : 'This will take around ', diag_est, ' s.'
781 : END IF
782 :
783 22 : CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
784 :
785 66 : ALLOCATE (Exc_ens(n_ov_joint))
786 :
787 22 : CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
788 :
789 22 : IF (diag_info /= 0) THEN
790 : CALL cp_abort(__LOCATION__, &
791 0 : "Diagonalization of A failed in TDA-BSE")
792 : END IF
793 :
794 22 : IF (nspins == 1) THEN
795 : CALL postprocess_bse(Exc_ens, fm_eigvec, bse_env, mo_coeff, &
796 14 : 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 8 : .TRUE., mo_coeff, bse_env, unit_nr)
800 : END IF
801 :
802 22 : CALL cp_fm_release(fm_eigvec)
803 22 : DEALLOCATE (Exc_ens)
804 :
805 22 : CALL timestop(handle)
806 :
807 66 : 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 12 : SUBROUTINE bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
825 12 : 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 12 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dip_ab_sp, fm_dip_ai_sp, fm_dip_ij_sp
849 48 : TYPE(cp_fm_type), DIMENSION(3) :: fm_trans_mom_joint
850 :
851 12 : CALL timeset(routineN, handle)
852 :
853 12 : nspins = SIZE(homo)
854 48 : ALLOCATE (n_ov_sp(nspins), offsets_sp(nspins))
855 12 : 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 12 : multiplet = "UKS"
860 12 : IF (flag_tda) THEN
861 8 : info_approximation = " -TDA- "
862 : ELSE
863 4 : info_approximation = "-ABBA-"
864 : END IF
865 :
866 12 : IF (unit_nr > 0) THEN
867 6 : 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 12 : 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 12 : 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 12 : CALL cp_fm_create(fm_trans_coeff, fm_eigvec_X%matrix_struct)
878 12 : CALL cp_fm_to_fm(fm_eigvec_X, fm_trans_coeff)
879 12 : 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 120 : ALLOCATE (fm_dip_ai_sp(3), fm_dip_ij_sp(3), fm_dip_ab_sp(3))
884 60 : ALLOCATE (oscill_str_joint(n_ov_joint), trans_mom_joint(3, 1, n_ov_joint))
885 36 : ALLOCATE (pol_res_joint(3, 3, n_ov_joint))
886 12 : trans_mom_joint(:, :, :) = 0.0_dp
887 12 : 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 12 : fm_trans_coeff%matrix_struct%context, 1, n_ov_joint)
890 48 : DO idir = 1, 3
891 36 : CALL cp_fm_create(fm_trans_mom_joint(idir), fm_struct_tmom)
892 48 : CALL cp_fm_set_all(fm_trans_mom_joint(idir), 0.0_dp)
893 : END DO
894 36 : 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 24 : ispin=isp)
899 24 : 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 24 : fm_trans_coeff%matrix_struct%context, n_ov_sp(isp), n_ov_joint)
902 24 : CALL cp_fm_create(fm_eigvec_sp, fm_struct_sp)
903 24 : 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 24 : 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 24 : fm_trans_coeff%matrix_struct%context, 1, n_ov_sp(isp))
908 96 : DO idir = 1, 3
909 72 : CALL cp_fm_create(fm_dip_reord_sp, fm_struct_dip_reord, name="bse_dip_reord")
910 72 : 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 72 : 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 72 : fm_dip_reord_sp, fm_eigvec_sp, 1.0_dp, fm_trans_mom_joint(idir))
915 72 : CALL cp_fm_release(fm_dip_reord_sp)
916 72 : CALL cp_fm_release(fm_dip_ai_sp(idir))
917 72 : CALL cp_fm_release(fm_dip_ij_sp(idir))
918 168 : CALL cp_fm_release(fm_dip_ab_sp(idir))
919 : END DO
920 24 : CALL cp_fm_release(fm_eigvec_sp)
921 24 : CALL cp_fm_struct_release(fm_struct_sp)
922 24 : NULLIFY (fm_struct_sp)
923 24 : CALL cp_fm_struct_release(fm_struct_dip_reord)
924 60 : NULLIFY (fm_struct_dip_reord)
925 : END DO
926 48 : DO idir = 1, 3
927 36 : CALL cp_fm_get_submatrix(fm_trans_mom_joint(idir), trans_mom_joint(idir, :, :))
928 48 : CALL cp_fm_release(fm_trans_mom_joint(idir))
929 : END DO
930 12 : CALL cp_fm_struct_release(fm_struct_tmom)
931 5692 : DO n = 1, n_ov_joint
932 22720 : DO idir = 1, 3
933 73840 : DO jdir = 1, 3
934 : pol_res_joint(idir, jdir, n) = 2.0_dp*Exc_ens(n)*trans_mom_joint(idir, 1, n) &
935 68160 : *trans_mom_joint(jdir, 1, n)
936 : END DO
937 : END DO
938 22732 : 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 36 : 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 12 : "implemented and have been skipped.")
947 12 : CALL cp_fm_release(fm_trans_coeff)
948 12 : DEALLOCATE (fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp)
949 12 : DEALLOCATE (n_ov_sp, offsets_sp, oscill_str_joint, trans_mom_joint, pol_res_joint)
950 :
951 12 : CALL timestop(handle)
952 :
953 48 : 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 38 : 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 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
986 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
987 : TYPE(cp_fm_type) :: fm_X_ia, fm_Y_ia
988 38 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
989 38 : fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
990 : TYPE(exciton_descr_type), ALLOCATABLE, &
991 38 : DIMENSION(:) :: exc_descr
992 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
993 :
994 38 : CALL timeset(routineN, handle)
995 :
996 : !Prepare variables for printing
997 38 : IF (bse_env%bse_spin_config == 0) THEN
998 38 : multiplet = "Singlet"
999 : ELSE
1000 0 : multiplet = "Triplet"
1001 : END IF
1002 38 : IF (.NOT. PRESENT(flag_TDA)) THEN
1003 0 : flag_TDA = .FALSE.
1004 : END IF
1005 38 : IF (flag_TDA) THEN
1006 16 : info_approximation = " -TDA- "
1007 : ELSE
1008 22 : info_approximation = "-ABBA-"
1009 : END IF
1010 :
1011 38 : n_moments_di = 3
1012 38 : 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 152 : ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
1016 152 : ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
1017 152 : ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
1018 38 : ALLOCATE (ref_point_multipole(3))
1019 152 : 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 38 : homo, virtual, fm_eigvec_X%matrix_struct%context)
1024 : ! Compute exciton descriptors from these multipoles
1025 38 : IF (bse_env%num_print_exc_descr > 0) THEN
1026 : ! Obtain quadrupoles in MO basis
1027 40 : ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
1028 40 : ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
1029 40 : ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
1030 : ! the preparer builds the AO quadrupoles whenever NUM_PRINT_EXC_DESCR is nonzero
1031 4 : 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 4 : homo, virtual, fm_eigvec_X%matrix_struct%context)
1035 : ! Iterate over excitation index outside of routine to make it compatible with tddft module
1036 424 : ALLOCATE (exc_descr(bse_env%num_print_exc_descr))
1037 104 : 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 100 : .FALSE., unit_nr, bse_env)
1040 100 : IF (.NOT. flag_TDA) THEN
1041 : CALL reshuffle_eigvec(fm_eigvec_Y, fm_Y_ia, homo, virtual, i_exc, &
1042 50 : .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 50 : 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 50 : i_exc, homo, virtual)
1054 : END IF
1055 100 : CALL cp_fm_release(fm_X_ia)
1056 104 : IF (.NOT. flag_TDA) THEN
1057 50 : CALL cp_fm_release(fm_Y_ia)
1058 : END IF
1059 : END DO
1060 : END IF
1061 :
1062 38 : 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 38 : fm_eigvec_Y)
1067 : END IF
1068 :
1069 : ! the iterative driver prints the header before its solver runs
1070 38 : IF (bse_env%bse_diag_method /= bse_iterdiag) THEN
1071 32 : 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 38 : 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 152 : 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 38 : info_approximation, bse_env, unit_nr)
1086 : ! Print exciton descriptors if keyword is invoked
1087 38 : IF (bse_env%num_print_exc_descr > 0) THEN
1088 4 : 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 4 : 'BSE|', particle_set)
1094 : END IF
1095 :
1096 : ! Compute and print excitation wavefunctions
1097 38 : IF (bse_env%do_nto_analysis) THEN
1098 4 : IF (unit_nr > 0) THEN
1099 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1100 : WRITE (unit_nr, '(T2,A4,T7,A47)') &
1101 2 : 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
1102 2 : 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 4 : bse_env, unit_nr)
1109 : END IF
1110 :
1111 152 : DO idir = 1, n_moments_di
1112 114 : CALL cp_fm_release(fm_dipole_ai_trunc(idir))
1113 114 : CALL cp_fm_release(fm_dipole_ij_trunc(idir))
1114 152 : CALL cp_fm_release(fm_dipole_ab_trunc(idir))
1115 : END DO
1116 38 : IF (bse_env%num_print_exc_descr > 0) THEN
1117 40 : DO idir = 1, n_moments_quad
1118 36 : CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
1119 36 : CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
1120 40 : CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
1121 : END DO
1122 4 : DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
1123 4 : DEALLOCATE (exc_descr)
1124 : END IF
1125 38 : DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
1126 38 : DEALLOCATE (ref_point_multipole)
1127 38 : IF (bse_env%bse_spin_config == 0) THEN
1128 38 : DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
1129 : END IF
1130 38 : IF (unit_nr > 0) THEN
1131 19 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1132 19 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1133 : END IF
1134 :
1135 38 : CALL timestop(handle)
1136 :
1137 76 : END SUBROUTINE postprocess_bse
1138 :
1139 : END MODULE bse_full_diag
|
