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