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