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