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: comp_eigvec_coeff_BSE,&
26 : fm_general_add_bse,&
27 : get_multipoles_mo,&
28 : reshuffle_eigvec
29 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
30 : cp_blacs_env_release,&
31 : cp_blacs_env_type
32 : USE cp_control_types, ONLY: dft_control_type,&
33 : tddfpt2_control_type
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
35 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
36 : cp_fm_power
37 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
38 : cp_fm_struct_release,&
39 : cp_fm_struct_type
40 : USE cp_fm_types, ONLY: cp_fm_create,&
41 : cp_fm_get_info,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE exstates_types, ONLY: excited_energy_type
47 : USE input_constants, ONLY: bse_screening_alpha,&
48 : bse_screening_rpa,&
49 : bse_singlet,&
50 : bse_triplet
51 : USE kinds, ONLY: dp
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE mp2_types, ONLY: mp2_type
54 : USE parallel_gemm_api, ONLY: parallel_gemm
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_environment_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag'
64 :
65 : PUBLIC :: create_A, diagonalize_A, create_B, create_hermitian_form_of_ABBA, &
66 : diagonalize_C
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
72 : !> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
73 : !> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
74 : !> α is a spin-dependent factor
75 : !> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
76 : !> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
77 : !> \param fm_mat_S_ia_bse ...
78 : !> \param fm_mat_S_bar_ij_bse ...
79 : !> \param fm_mat_S_ab_bse ...
80 : !> \param fm_A ...
81 : !> \param Eigenval ...
82 : !> \param unit_nr ...
83 : !> \param homo ...
84 : !> \param virtual ...
85 : !> \param dimen_RI ...
86 : !> \param mp2_env ...
87 : !> \param para_env ...
88 : !> \param qs_env ...
89 : ! **************************************************************************************************
90 170 : SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
91 34 : fm_A, Eigenval, unit_nr, &
92 : homo, virtual, dimen_RI, mp2_env, &
93 : para_env, qs_env)
94 :
95 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, &
96 : fm_mat_S_ab_bse
97 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
98 : REAL(KIND=dp), DIMENSION(:) :: Eigenval
99 : INTEGER, INTENT(IN) :: unit_nr, homo, virtual, dimen_RI
100 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
101 : TYPE(mp_para_env_type), INTENT(INOUT) :: para_env
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 :
104 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_A'
105 :
106 : INTEGER :: a_virt_row, handle, i_occ_row, &
107 : i_row_global, ii, j_col_global, jj, &
108 : ncol_local_A, nrow_local_A, sizeeigen
109 : INTEGER, DIMENSION(4) :: reordering
110 34 : INTEGER, DIMENSION(:), POINTER :: col_indices_A, row_indices_A
111 : REAL(KIND=dp) :: alpha, alpha_screening, eigen_diff
112 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
113 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_A, fm_struct_W
114 : TYPE(cp_fm_type) :: fm_A_copy, fm_W
115 : TYPE(dft_control_type), POINTER :: dft_control
116 : TYPE(excited_energy_type), POINTER :: ex_env
117 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
118 :
119 34 : CALL timeset(routineN, handle)
120 :
121 34 : NULLIFY (dft_control, tddfpt_control)
122 34 : CALL get_qs_env(qs_env, dft_control=dft_control)
123 34 : tddfpt_control => dft_control%tddfpt2_control
124 :
125 34 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
126 3 : WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
127 : END IF
128 :
129 : !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
130 68 : SELECT CASE (mp2_env%bse%bse_spin_config)
131 : CASE (bse_singlet)
132 34 : alpha = 2.0_dp
133 : CASE (bse_triplet)
134 34 : alpha = 0.0_dp
135 : END SELECT
136 :
137 34 : IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
138 2 : alpha_screening = mp2_env%bse%screening_factor
139 : ELSE
140 32 : alpha_screening = 1.0_dp
141 : END IF
142 :
143 : ! create the blacs env for ij matrices (NOT fm_mat_S_ia_bse%matrix_struct related parallel_gemms!)
144 34 : NULLIFY (blacs_env)
145 34 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
146 :
147 : ! We have to use the same blacs_env for A as for the matrices fm_mat_S_ia_bse from RPA
148 : ! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
149 : ! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab
150 : ! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb
151 : ! We use the A matrix already from the start instead of v
152 : CALL cp_fm_struct_create(fm_struct_A, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
153 34 : ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
154 34 : CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb")
155 34 : CALL cp_fm_set_all(fm_A, 0.0_dp)
156 34 : IF (tddfpt_control%do_bse_w_only) THEN
157 2 : CALL cp_fm_create(fm_A_copy, fm_struct_A, name="fm_A_iajb")
158 2 : CALL cp_fm_set_all(fm_A_copy, 0.0_dp)
159 : END IF
160 :
161 : CALL cp_fm_struct_create(fm_struct_W, context=fm_mat_S_ab_bse%matrix_struct%context, nrow_global=homo**2, &
162 34 : ncol_global=virtual**2, para_env=fm_mat_S_ab_bse%matrix_struct%para_env)
163 34 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab")
164 34 : CALL cp_fm_set_all(fm_W, 0.0_dp)
165 :
166 : ! Create A matrix from GW Energies, v_ia,jb and W_ij,ab (different blacs_env!)
167 : ! v_ia,jb, which is directly initialized in A (with a factor of alpha)
168 : ! v_ia,jb = \sum_P B^P_ia B^P_jb
169 34 : IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN
170 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
171 : matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
172 30 : matrix_c=fm_A)
173 : END IF
174 :
175 34 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
176 3 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
177 : END IF
178 :
179 : ! If infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
180 34 : IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
181 : !W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab
182 : CALL parallel_gemm(transa="T", transb="N", m=homo**2, n=virtual**2, k=dimen_RI, alpha=alpha_screening, &
183 : matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, &
184 32 : matrix_c=fm_W)
185 : END IF
186 :
187 34 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
188 3 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
189 : END IF
190 :
191 : ! We start by moving data from local parts of W_ij,ab to the full matrix A_ia,jb using buffers
192 : CALL cp_fm_get_info(matrix=fm_A, &
193 : nrow_local=nrow_local_A, &
194 : ncol_local=ncol_local_A, &
195 : row_indices=row_indices_A, &
196 34 : col_indices=col_indices_A)
197 : ! Writing -1.0_dp * W_ij,ab to A_ia,jb, i.e. beta = -1.0_dp,
198 : ! W_ij,ab: nrow_secidx_in = homo, ncol_secidx_in = virtual
199 : ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
200 :
201 : ! If infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
202 34 : IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
203 32 : reordering = [1, 3, 2, 4]
204 : CALL fm_general_add_bse(fm_A, fm_W, -1.0_dp, homo, virtual, &
205 32 : virtual, virtual, unit_nr, reordering, mp2_env)
206 32 : IF (tddfpt_control%do_bse_w_only) THEN
207 : CALL fm_general_add_bse(fm_A_copy, fm_W, -1.0_dp, homo, virtual, &
208 2 : virtual, virtual, unit_nr, reordering, mp2_env)
209 : END IF
210 : END IF
211 : !full matrix W is not needed anymore, release it to save memory
212 34 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
213 : tddfpt_control%do_bse_gw_only) THEN
214 4 : NULLIFY (ex_env)
215 4 : CALL get_qs_env(qs_env, exstate_env=ex_env)
216 4 : IF (.NOT. tddfpt_control%do_bse_gw_only) THEN
217 12 : ALLOCATE (ex_env%bse_w_matrix_MO(1, 1)) ! for now only closed-shell
218 12 : ALLOCATE (ex_env%bse_a_matrix_MO(1, 1)) ! for now only closed-shell
219 4 : CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_W)
220 4 : CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_A)
221 4 : CALL cp_fm_to_fm(fm_W, ex_env%bse_w_matrix_MO(1, 1))
222 4 : IF (tddfpt_control%do_bse_w_only) THEN
223 2 : CALL cp_fm_to_fm(fm_A_copy, ex_env%bse_a_matrix_MO(1, 1))
224 : ELSE
225 2 : CALL cp_fm_to_fm(fm_A, ex_env%bse_a_matrix_MO(1, 1))
226 : END IF
227 : END IF
228 : END IF
229 34 : CALL cp_fm_release(fm_W)
230 34 : IF (tddfpt_control%do_bse_w_only) CALL cp_fm_release(fm_A_copy)
231 :
232 : !Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb
233 34 : IF (.NOT. tddfpt_control%do_bse) THEN
234 788 : DO ii = 1, nrow_local_A
235 :
236 756 : i_row_global = row_indices_A(ii)
237 :
238 38884 : DO jj = 1, ncol_local_A
239 :
240 38096 : j_col_global = col_indices_A(jj)
241 :
242 38852 : IF (i_row_global == j_col_global) THEN
243 756 : i_occ_row = (i_row_global - 1)/virtual + 1
244 756 : a_virt_row = MOD(i_row_global - 1, virtual) + 1
245 756 : eigen_diff = Eigenval(a_virt_row + homo) - Eigenval(i_occ_row)
246 756 : fm_A%local_data(ii, jj) = fm_A%local_data(ii, jj) + eigen_diff
247 :
248 : END IF
249 : END DO
250 : END DO
251 : END IF
252 :
253 34 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. tddfpt_control%do_bse_gw_only) THEN
254 4 : sizeeigen = SIZE(Eigenval)
255 12 : ALLOCATE (ex_env%gw_eigen(sizeeigen)) ! for now only closed-shell
256 126 : ex_env%gw_eigen(:) = Eigenval(:)
257 : END IF
258 :
259 34 : CALL cp_fm_struct_release(fm_struct_A)
260 34 : CALL cp_fm_struct_release(fm_struct_W)
261 :
262 34 : CALL cp_blacs_env_release(blacs_env)
263 :
264 34 : CALL timestop(handle)
265 :
266 34 : END SUBROUTINE create_A
267 :
268 : ! **************************************************************************************************
269 : !> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
270 : !> B_ia,jb = α * v_ia,jb - W_ib,aj
271 : !> α is a spin-dependent factor
272 : !> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
273 : !> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
274 : !> \param fm_mat_S_ia_bse ...
275 : !> \param fm_mat_S_bar_ia_bse ...
276 : !> \param fm_B ...
277 : !> \param homo ...
278 : !> \param virtual ...
279 : !> \param dimen_RI ...
280 : !> \param unit_nr ...
281 : !> \param mp2_env ...
282 : ! **************************************************************************************************
283 54 : SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
284 : homo, virtual, dimen_RI, unit_nr, mp2_env)
285 :
286 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse
287 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_B
288 : INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr
289 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
290 :
291 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_B'
292 :
293 : INTEGER :: handle
294 : INTEGER, DIMENSION(4) :: reordering
295 : REAL(KIND=dp) :: alpha, alpha_screening
296 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_v
297 : TYPE(cp_fm_type) :: fm_W
298 :
299 18 : CALL timeset(routineN, handle)
300 :
301 18 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
302 2 : WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
303 : END IF
304 :
305 : ! Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
306 36 : SELECT CASE (mp2_env%bse%bse_spin_config)
307 : CASE (bse_singlet)
308 18 : alpha = 2.0_dp
309 : CASE (bse_triplet)
310 18 : alpha = 0.0_dp
311 : END SELECT
312 :
313 18 : IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
314 2 : alpha_screening = mp2_env%bse%screening_factor
315 : ELSE
316 16 : alpha_screening = 1.0_dp
317 : END IF
318 :
319 : CALL cp_fm_struct_create(fm_struct_v, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
320 18 : ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
321 18 : CALL cp_fm_create(fm_B, fm_struct_v, name="fm_B_iajb")
322 18 : CALL cp_fm_set_all(fm_B, 0.0_dp)
323 :
324 18 : CALL cp_fm_create(fm_W, fm_struct_v, name="fm_W_ibaj")
325 18 : CALL cp_fm_set_all(fm_W, 0.0_dp)
326 :
327 18 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
328 2 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
329 : END IF
330 : ! v_ia,jb = \sum_P B^P_ia B^P_jb
331 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
332 : matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
333 18 : matrix_c=fm_B)
334 :
335 : ! If infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
336 18 : IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
337 : ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj
338 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha_screening, &
339 : matrix_a=fm_mat_S_bar_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
340 16 : matrix_c=fm_W)
341 :
342 : ! from W_ib,ja to A_ia,jb (formally: W_ib,aj, but our internal indexorder is different)
343 : ! Writing -1.0_dp * W_ib,ja to A_ia,jb, i.e. beta = -1.0_dp,
344 : ! W_ib,ja: nrow_secidx_in = virtual, ncol_secidx_in = virtual
345 : ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
346 16 : reordering = [1, 4, 3, 2]
347 : CALL fm_general_add_bse(fm_B, fm_W, -1.0_dp, virtual, virtual, &
348 16 : virtual, virtual, unit_nr, reordering, mp2_env)
349 : END IF
350 :
351 18 : CALL cp_fm_release(fm_W)
352 18 : CALL cp_fm_struct_release(fm_struct_v)
353 18 : CALL timestop(handle)
354 :
355 18 : END SUBROUTINE create_B
356 :
357 : ! **************************************************************************************************
358 : !> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
359 : !> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
360 : !> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
361 : !> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
362 : !> \param fm_A ...
363 : !> \param fm_B ...
364 : !> \param fm_C ...
365 : !> \param fm_sqrt_A_minus_B ...
366 : !> \param fm_inv_sqrt_A_minus_B ...
367 : !> \param homo ...
368 : !> \param virtual ...
369 : !> \param unit_nr ...
370 : !> \param mp2_env ...
371 : !> \param diag_est ...
372 : ! **************************************************************************************************
373 144 : SUBROUTINE create_hermitian_form_of_ABBA(fm_A, fm_B, fm_C, &
374 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
375 : homo, virtual, unit_nr, mp2_env, diag_est)
376 :
377 : TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B
378 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C, fm_sqrt_A_minus_B, &
379 : fm_inv_sqrt_A_minus_B
380 : INTEGER, INTENT(IN) :: homo, virtual, unit_nr
381 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
382 : REAL(KIND=dp), INTENT(IN) :: diag_est
383 :
384 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_hermitian_form_of_ABBA'
385 :
386 : INTEGER :: dim_mat, handle, n_dependent
387 : REAL(KIND=dp), DIMENSION(2) :: eigvals_AB_diff
388 : TYPE(cp_fm_type) :: fm_A_minus_B, fm_A_plus_B, fm_dummy, &
389 : fm_work_product
390 :
391 18 : CALL timeset(routineN, handle)
392 :
393 18 : IF (unit_nr > 0) THEN
394 9 : WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
395 18 : ' with size of A. This will take around ', diag_est, " s."
396 : END IF
397 :
398 : ! Create work matrices, which will hold A+B and A-B and their powers
399 : ! C is created afterwards to save memory
400 : ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
401 : ! \_______/ \___/ \______/
402 : ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
403 : ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
404 : ! Intermediate work matrices:
405 : ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
406 : ! fm_A_minus_B: (A-B) EQ.III
407 : ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
408 18 : CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct)
409 18 : CALL cp_fm_to_fm(fm_A, fm_A_plus_B)
410 18 : CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct)
411 18 : CALL cp_fm_to_fm(fm_A, fm_A_minus_B)
412 18 : CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct)
413 18 : CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp)
414 18 : CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct)
415 18 : CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp)
416 :
417 18 : CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
418 :
419 18 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
420 2 : WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
421 : END IF
422 :
423 : ! Add/Substract B (cf. EQs. Ib and III)
424 18 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B)
425 18 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_minus_B, -1.0_dp, fm_B)
426 :
427 : ! cp_fm_power will overwrite matrix, therefore we create copies
428 18 : CALL cp_fm_to_fm(fm_A_minus_B, fm_inv_sqrt_A_minus_B)
429 :
430 : ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
431 : ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
432 :
433 : ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
434 18 : CALL cp_fm_create(fm_dummy, fm_A%matrix_struct)
435 : ! Create (A-B)^-0.5 (cf. EQ.II)
436 18 : CALL cp_fm_power(fm_inv_sqrt_A_minus_B, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_AB_diff)
437 18 : CALL cp_fm_release(fm_dummy)
438 : ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
439 : ! In this case, the procedure for hermitian form of ABBA is not applicable
440 18 : IF (eigvals_AB_diff(1) < 0) THEN
441 : CALL cp_abort(__LOCATION__, &
442 : "Matrix (A-B) is not positive definite. "// &
443 0 : "Hermitian diagonalization of full ABBA matrix is ill-defined.")
444 : END IF
445 :
446 : ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
447 : ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
448 : ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
449 18 : dim_mat = homo*virtual
450 : 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, &
451 18 : fm_sqrt_A_minus_B)
452 :
453 : ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
454 : 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, &
455 18 : fm_work_product)
456 :
457 : ! Release to save memory
458 18 : CALL cp_fm_release(fm_A_plus_B)
459 18 : CALL cp_fm_release(fm_A_minus_B)
460 :
461 : ! Now create full
462 18 : CALL cp_fm_create(fm_C, fm_A%matrix_struct)
463 18 : CALL cp_fm_set_all(fm_C, 0.0_dp)
464 : ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
465 : 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, &
466 18 : fm_C)
467 18 : CALL cp_fm_release(fm_work_product)
468 :
469 18 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
470 2 : WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
471 : END IF
472 :
473 18 : CALL timestop(handle)
474 18 : END SUBROUTINE create_hermitian_form_of_ABBA
475 :
476 : ! **************************************************************************************************
477 : !> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
478 : !> Here, the eigenvectors Z^n relate to X^n via
479 : !> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
480 : !> \param fm_C ...
481 : !> \param homo ...
482 : !> \param virtual ...
483 : !> \param homo_irred ...
484 : !> \param fm_sqrt_A_minus_B ...
485 : !> \param fm_inv_sqrt_A_minus_B ...
486 : !> \param unit_nr ...
487 : !> \param diag_est ...
488 : !> \param mp2_env ...
489 : !> \param qs_env ...
490 : !> \param mo_coeff ...
491 : ! **************************************************************************************************
492 18 : SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
493 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
494 18 : unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
495 :
496 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C
497 : INTEGER, INTENT(IN) :: homo, virtual, homo_irred
498 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B
499 : INTEGER, INTENT(IN) :: unit_nr
500 : REAL(KIND=dp), INTENT(IN) :: diag_est
501 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
502 : TYPE(qs_environment_type), POINTER :: qs_env
503 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
504 :
505 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_C'
506 :
507 : INTEGER :: diag_info, handle
508 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
509 : TYPE(cp_fm_type) :: fm_eigvec_X, fm_eigvec_Y, fm_eigvec_Z, &
510 : fm_mat_eigvec_transform_diff, &
511 : fm_mat_eigvec_transform_sum
512 :
513 18 : CALL timeset(routineN, handle)
514 :
515 18 : IF (unit_nr > 0) THEN
516 9 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
517 18 : 'This will take around ', diag_est, ' s.'
518 : END IF
519 :
520 : !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
521 : !Now: Diagonalize it
522 18 : CALL cp_fm_create(fm_eigvec_Z, fm_C%matrix_struct)
523 :
524 54 : ALLOCATE (Exc_ens(homo*virtual))
525 :
526 18 : CALL choose_eigv_solver(fm_C, fm_eigvec_Z, Exc_ens, diag_info)
527 :
528 18 : IF (diag_info /= 0) THEN
529 : CALL cp_abort(__LOCATION__, &
530 0 : "Diagonalization of C=(A-B)^0.5 (A+B) (A-B)^0.5 failed in BSE")
531 : END IF
532 :
533 : ! C could have negative eigenvalues, since we do not explicitly check A+B
534 : ! for positive definiteness (would make another O(N^6) Diagon. necessary)
535 : ! Instead, we include a check here
536 18 : IF (Exc_ens(1) < 0) THEN
537 0 : IF (unit_nr > 0) THEN
538 : CALL cp_abort(__LOCATION__, &
539 : "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
540 0 : "(A+B) is not positive definite.")
541 : END IF
542 : END IF
543 882 : Exc_ens = SQRT(Exc_ens)
544 :
545 : ! Prepare eigenvector for interpretation of singleparticle transitions
546 : ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
547 : ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
548 :
549 : ! Following Furche, we basically use Eqs. (A10): First, we multiply
550 : ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
551 : ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
552 :
553 : ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
554 18 : CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_C%matrix_struct)
555 18 : CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
556 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
557 : matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
558 18 : matrix_c=fm_mat_eigvec_transform_sum)
559 18 : CALL cp_fm_release(fm_sqrt_A_minus_B)
560 : ! This normalizes the eigenvectors
561 18 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_sum, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
562 :
563 : ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
564 18 : CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_C%matrix_struct)
565 18 : CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
566 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
567 : matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
568 18 : matrix_c=fm_mat_eigvec_transform_diff)
569 18 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
570 18 : CALL cp_fm_release(fm_eigvec_Z)
571 :
572 : ! This normalizes the eigenvectors
573 18 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_diff, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
574 :
575 : ! Now, we add the two equations to obtain X_n
576 : ! Add overwrites the first argument, therefore we copy it beforehand
577 18 : CALL cp_fm_create(fm_eigvec_X, fm_C%matrix_struct)
578 18 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_X)
579 18 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_X, 1.0_dp, fm_mat_eigvec_transform_diff)
580 :
581 : ! Now, we subtract the two equations to obtain Y_n
582 : ! Add overwrites the first argument, therefore we copy it beforehand
583 18 : CALL cp_fm_create(fm_eigvec_Y, fm_C%matrix_struct)
584 18 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_Y)
585 18 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_Y, -1.0_dp, fm_mat_eigvec_transform_diff)
586 :
587 : !Cleanup
588 18 : CALL cp_fm_release(fm_mat_eigvec_transform_diff)
589 18 : CALL cp_fm_release(fm_mat_eigvec_transform_sum)
590 :
591 : CALL postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
592 : homo, virtual, homo_irred, unit_nr, &
593 18 : .FALSE., fm_eigvec_Y)
594 :
595 18 : DEALLOCATE (Exc_ens)
596 18 : CALL cp_fm_release(fm_eigvec_X)
597 18 : CALL cp_fm_release(fm_eigvec_Y)
598 :
599 18 : CALL timestop(handle)
600 :
601 126 : END SUBROUTINE diagonalize_C
602 :
603 : ! **************************************************************************************************
604 : !> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
605 : !> \param fm_A ...
606 : !> \param homo ...
607 : !> \param virtual ...
608 : !> \param homo_irred ...
609 : !> \param unit_nr ...
610 : !> \param diag_est ...
611 : !> \param mp2_env ...
612 : !> \param qs_env ...
613 : !> \param mo_coeff ...
614 : ! **************************************************************************************************
615 14 : SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
616 14 : unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
617 :
618 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
619 : INTEGER, INTENT(IN) :: homo, virtual, homo_irred, unit_nr
620 : REAL(KIND=dp), INTENT(IN) :: diag_est
621 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
622 : TYPE(qs_environment_type), POINTER :: qs_env
623 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
624 :
625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A'
626 :
627 : INTEGER :: diag_info, handle
628 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
629 : TYPE(cp_fm_type) :: fm_eigvec
630 :
631 14 : CALL timeset(routineN, handle)
632 :
633 : !Continue with formatting of subroutine create_A
634 14 : IF (unit_nr > 0) THEN
635 7 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
636 14 : 'This will take around ', diag_est, ' s.'
637 : END IF
638 :
639 : !We have now the full matrix A_iajb, distributed over all ranks
640 : !Now: Diagonalize it
641 14 : CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
642 :
643 42 : ALLOCATE (Exc_ens(homo*virtual))
644 :
645 14 : CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
646 :
647 14 : IF (diag_info /= 0) THEN
648 : CALL cp_abort(__LOCATION__, &
649 0 : "Diagonalization of A failed in TDA-BSE")
650 : END IF
651 :
652 : CALL postprocess_bse(Exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
653 14 : homo, virtual, homo_irred, unit_nr, .TRUE.)
654 :
655 14 : CALL cp_fm_release(fm_eigvec)
656 14 : DEALLOCATE (Exc_ens)
657 :
658 14 : CALL timestop(handle)
659 :
660 42 : END SUBROUTINE diagonalize_A
661 :
662 : ! **************************************************************************************************
663 : !> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
664 : !> \param Exc_ens ...
665 : !> \param fm_eigvec_X ...
666 : !> \param mp2_env ...
667 : !> \param qs_env ...
668 : !> \param mo_coeff ...
669 : !> \param homo ...
670 : !> \param virtual ...
671 : !> \param homo_irred ...
672 : !> \param unit_nr ...
673 : !> \param flag_TDA ...
674 : !> \param fm_eigvec_Y ...
675 : ! **************************************************************************************************
676 32 : SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
677 : homo, virtual, homo_irred, unit_nr, &
678 : flag_TDA, fm_eigvec_Y)
679 :
680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
681 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
682 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
683 : TYPE(qs_environment_type), POINTER :: qs_env
684 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
685 : INTEGER :: homo, virtual, homo_irred, unit_nr
686 : LOGICAL, OPTIONAL :: flag_TDA
687 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
688 :
689 : CHARACTER(LEN=*), PARAMETER :: routineN = 'postprocess_bse'
690 :
691 : CHARACTER(LEN=10) :: info_approximation, multiplet
692 : INTEGER :: handle, i_exc, idir, n_moments_di, &
693 : n_moments_quad
694 : REAL(KIND=dp) :: alpha
695 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
696 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
697 : TYPE(cp_fm_type) :: fm_X_ia, fm_Y_ia
698 32 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
699 32 : fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
700 : TYPE(exciton_descr_type), ALLOCATABLE, &
701 32 : DIMENSION(:) :: exc_descr
702 :
703 32 : CALL timeset(routineN, handle)
704 :
705 : !Prepare variables for printing
706 32 : IF (mp2_env%bse%bse_spin_config == 0) THEN
707 32 : multiplet = "Singlet"
708 32 : alpha = 2.0_dp
709 : ELSE
710 0 : multiplet = "Triplet"
711 0 : alpha = 0.0_dp
712 : END IF
713 32 : IF (.NOT. PRESENT(flag_TDA)) THEN
714 0 : flag_TDA = .FALSE.
715 : END IF
716 32 : IF (flag_TDA) THEN
717 14 : info_approximation = " -TDA- "
718 : ELSE
719 18 : info_approximation = "-ABBA-"
720 : END IF
721 :
722 32 : n_moments_di = 3
723 32 : n_moments_quad = 9
724 : ! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
725 : ! Need dipoles also for spatial expectation values, which are well-defined also for triplets
726 128 : ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
727 128 : ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
728 128 : ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
729 32 : ALLOCATE (ref_point_multipole(3))
730 : ! Obtain dipoles in MO basis
731 : CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
732 : qs_env, mo_coeff, ref_point_multipole, 1, &
733 32 : homo, virtual, fm_eigvec_X%matrix_struct%context)
734 : ! Compute exciton descriptors from these multipoles
735 32 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
736 : ! Obtain quadrupoles in MO basis
737 40 : ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
738 40 : ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
739 40 : ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
740 : CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
741 : qs_env, mo_coeff, ref_point_multipole, 2, &
742 4 : homo, virtual, fm_eigvec_X%matrix_struct%context)
743 : ! Iterate over excitation index outside of routine to make it compatible with tddft module
744 424 : ALLOCATE (exc_descr(mp2_env%bse%num_print_exc_descr))
745 104 : DO i_exc = 1, mp2_env%bse%num_print_exc_descr
746 : CALL reshuffle_eigvec(fm_eigvec_X, fm_X_ia, homo, virtual, i_exc, &
747 100 : .FALSE., unit_nr, mp2_env)
748 100 : IF (.NOT. flag_TDA) THEN
749 : CALL reshuffle_eigvec(fm_eigvec_Y, fm_Y_ia, homo, virtual, i_exc, &
750 50 : .FALSE., unit_nr, mp2_env)
751 :
752 : CALL get_exciton_descriptors(exc_descr, fm_X_ia, &
753 : fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
754 : fm_quadpole_ai_trunc, &
755 : i_exc, homo, virtual, &
756 50 : fm_Y_ia)
757 : ELSE
758 : CALL get_exciton_descriptors(exc_descr, fm_X_ia, &
759 : fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
760 : fm_quadpole_ai_trunc, &
761 50 : i_exc, homo, virtual)
762 : END IF
763 100 : CALL cp_fm_release(fm_X_ia)
764 104 : IF (.NOT. flag_TDA) THEN
765 50 : CALL cp_fm_release(fm_Y_ia)
766 : END IF
767 : END DO
768 : END IF
769 :
770 32 : IF (mp2_env%bse%bse_spin_config == 0) THEN
771 : CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
772 : trans_mom_bse, oscill_str, polarizability_residues, &
773 : mp2_env, homo, virtual, unit_nr, &
774 32 : fm_eigvec_Y)
775 : END IF
776 :
777 : ! Prints basic definitions used in BSE calculation
778 : CALL print_output_header(homo, virtual, homo_irred, flag_TDA, &
779 32 : multiplet, alpha, mp2_env, unit_nr)
780 :
781 : ! Prints excitation energies up to user-specified number
782 : CALL print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
783 32 : info_approximation, mp2_env, unit_nr)
784 :
785 : ! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
786 : CALL print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
787 32 : info_approximation, mp2_env, unit_nr, fm_eigvec_Y)
788 :
789 : ! Prints optical properties, if state is a singlet
790 : CALL print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
791 : homo, virtual, homo_irred, flag_TDA, &
792 32 : info_approximation, mp2_env, unit_nr)
793 : ! Print exciton descriptors if keyword is invoked
794 32 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
795 : CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
796 : mp2_env%bse%num_print_exc_descr, mp2_env%bse%bse_debug_print, &
797 : mp2_env%bse%print_directional_exc_descr, &
798 4 : 'BSE|', qs_env)
799 : END IF
800 :
801 : ! Compute and print excitation wavefunctions
802 32 : IF (mp2_env%bse%do_nto_analysis) THEN
803 4 : IF (unit_nr > 0) THEN
804 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
805 : WRITE (unit_nr, '(T2,A4,T7,A47)') &
806 2 : 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
807 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
808 : END IF
809 : CALL calculate_NTOs(fm_eigvec_X, fm_eigvec_Y, &
810 : mo_coeff, homo, virtual, &
811 : info_approximation, &
812 : oscill_str, &
813 4 : qs_env, unit_nr, mp2_env)
814 : END IF
815 :
816 128 : DO idir = 1, n_moments_di
817 96 : CALL cp_fm_release(fm_dipole_ai_trunc(idir))
818 96 : CALL cp_fm_release(fm_dipole_ij_trunc(idir))
819 128 : CALL cp_fm_release(fm_dipole_ab_trunc(idir))
820 : END DO
821 32 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
822 40 : DO idir = 1, n_moments_quad
823 36 : CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
824 36 : CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
825 40 : CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
826 : END DO
827 4 : DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
828 4 : DEALLOCATE (exc_descr)
829 : END IF
830 32 : DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
831 32 : DEALLOCATE (ref_point_multipole)
832 32 : IF (mp2_env%bse%bse_spin_config == 0) THEN
833 32 : DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
834 : END IF
835 32 : IF (unit_nr > 0) THEN
836 16 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
837 16 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
838 : END IF
839 :
840 32 : CALL timestop(handle)
841 :
842 64 : END SUBROUTINE postprocess_bse
843 :
844 : END MODULE bse_full_diag
|
