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 26 : CPASSERT(diag_info == 0)
479 : ! C could have negative eigenvalues, since we do not explicitly check A+B
480 : ! for positive definiteness (would make another O(N^6) Diagon. necessary)
481 : ! Instead, we include a check here
482 26 : IF (Exc_ens(1) < 0) THEN
483 0 : IF (unit_nr > 0) THEN
484 : CALL cp_abort(__LOCATION__, &
485 : "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
486 0 : "(A+B) is not positive definite.")
487 : END IF
488 : END IF
489 1274 : Exc_ens = SQRT(Exc_ens)
490 :
491 : ! Prepare eigenvector for interpretation of singleparticle transitions
492 : ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
493 : ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
494 :
495 : ! Following Furche, we basically use Eqs. (A10): First, we multiply
496 : ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
497 : ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
498 :
499 : ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
500 26 : CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_C%matrix_struct)
501 26 : CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
502 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
503 : matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
504 26 : matrix_c=fm_mat_eigvec_transform_sum)
505 26 : CALL cp_fm_release(fm_sqrt_A_minus_B)
506 : ! This normalizes the eigenvectors
507 26 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_sum, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
508 :
509 : ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
510 26 : CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_C%matrix_struct)
511 26 : CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
512 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
513 : matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
514 26 : matrix_c=fm_mat_eigvec_transform_diff)
515 26 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
516 26 : CALL cp_fm_release(fm_eigvec_Z)
517 :
518 : ! This normalizes the eigenvectors
519 26 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_diff, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
520 :
521 : ! Now, we add the two equations to obtain X_n
522 : ! Add overwrites the first argument, therefore we copy it beforehand
523 26 : CALL cp_fm_create(fm_eigvec_X, fm_C%matrix_struct)
524 26 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_X)
525 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_X, 1.0_dp, fm_mat_eigvec_transform_diff)
526 :
527 : ! Now, we subtract the two equations to obtain Y_n
528 : ! Add overwrites the first argument, therefore we copy it beforehand
529 26 : CALL cp_fm_create(fm_eigvec_Y, fm_C%matrix_struct)
530 26 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_Y)
531 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_Y, -1.0_dp, fm_mat_eigvec_transform_diff)
532 :
533 : !Cleanup
534 26 : CALL cp_fm_release(fm_mat_eigvec_transform_diff)
535 26 : CALL cp_fm_release(fm_mat_eigvec_transform_sum)
536 :
537 : CALL postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
538 : homo, virtual, homo_irred, unit_nr, &
539 26 : .FALSE., fm_eigvec_Y)
540 :
541 26 : DEALLOCATE (Exc_ens)
542 26 : CALL cp_fm_release(fm_eigvec_X)
543 26 : CALL cp_fm_release(fm_eigvec_Y)
544 :
545 26 : CALL timestop(handle)
546 :
547 182 : END SUBROUTINE diagonalize_C
548 :
549 : ! **************************************************************************************************
550 : !> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
551 : !> \param fm_A ...
552 : !> \param homo ...
553 : !> \param virtual ...
554 : !> \param homo_irred ...
555 : !> \param unit_nr ...
556 : !> \param diag_est ...
557 : !> \param mp2_env ...
558 : !> \param qs_env ...
559 : !> \param mo_coeff ...
560 : ! **************************************************************************************************
561 16 : SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
562 16 : unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
563 :
564 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
565 : INTEGER, INTENT(IN) :: homo, virtual, homo_irred, unit_nr
566 : REAL(KIND=dp), INTENT(IN) :: diag_est
567 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
568 : TYPE(qs_environment_type), POINTER :: qs_env
569 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
570 :
571 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A'
572 :
573 : INTEGER :: diag_info, handle
574 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
575 : TYPE(cp_fm_type) :: fm_eigvec
576 :
577 16 : CALL timeset(routineN, handle)
578 :
579 : !Continue with formatting of subroutine create_A
580 16 : IF (unit_nr > 0) THEN
581 8 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
582 16 : 'This will take around ', diag_est, ' s.'
583 : END IF
584 :
585 : !We have now the full matrix A_iajb, distributed over all ranks
586 : !Now: Diagonalize it
587 16 : CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
588 :
589 48 : ALLOCATE (Exc_ens(homo*virtual))
590 :
591 16 : CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
592 16 : CPASSERT(diag_info == 0)
593 : CALL postprocess_bse(Exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
594 16 : homo, virtual, homo_irred, unit_nr, .TRUE.)
595 :
596 16 : CALL cp_fm_release(fm_eigvec)
597 16 : DEALLOCATE (Exc_ens)
598 :
599 16 : CALL timestop(handle)
600 :
601 48 : END SUBROUTINE diagonalize_A
602 :
603 : ! **************************************************************************************************
604 : !> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
605 : !> \param Exc_ens ...
606 : !> \param fm_eigvec_X ...
607 : !> \param mp2_env ...
608 : !> \param qs_env ...
609 : !> \param mo_coeff ...
610 : !> \param homo ...
611 : !> \param virtual ...
612 : !> \param homo_irred ...
613 : !> \param unit_nr ...
614 : !> \param flag_TDA ...
615 : !> \param fm_eigvec_Y ...
616 : ! **************************************************************************************************
617 42 : SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
618 : homo, virtual, homo_irred, unit_nr, &
619 : flag_TDA, fm_eigvec_Y)
620 :
621 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
622 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
623 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
624 : TYPE(qs_environment_type), POINTER :: qs_env
625 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
626 : INTEGER :: homo, virtual, homo_irred, unit_nr
627 : LOGICAL, OPTIONAL :: flag_TDA
628 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
629 :
630 : CHARACTER(LEN=*), PARAMETER :: routineN = 'postprocess_bse'
631 :
632 : CHARACTER(LEN=10) :: info_approximation, multiplet
633 : INTEGER :: handle, i_exc, idir, n_moments_di, &
634 : n_moments_quad
635 : REAL(KIND=dp) :: alpha
636 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
637 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
638 : TYPE(cp_fm_type) :: fm_X_ia, fm_Y_ia
639 42 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
640 42 : fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
641 : TYPE(exciton_descr_type), ALLOCATABLE, &
642 42 : DIMENSION(:) :: exc_descr
643 :
644 42 : CALL timeset(routineN, handle)
645 :
646 : !Prepare variables for printing
647 42 : IF (mp2_env%bse%bse_spin_config == 0) THEN
648 42 : multiplet = "Singlet"
649 42 : alpha = 2.0_dp
650 : ELSE
651 0 : multiplet = "Triplet"
652 0 : alpha = 0.0_dp
653 : END IF
654 42 : IF (.NOT. PRESENT(flag_TDA)) THEN
655 0 : flag_TDA = .FALSE.
656 : END IF
657 42 : IF (flag_TDA) THEN
658 16 : info_approximation = " -TDA- "
659 : ELSE
660 26 : info_approximation = "-ABBA-"
661 : END IF
662 :
663 42 : n_moments_di = 3
664 42 : n_moments_quad = 9
665 : ! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
666 : ! Need dipoles also for spatial expectation values, which are well-defined also for triplets
667 168 : ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
668 168 : ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
669 168 : ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
670 42 : ALLOCATE (ref_point_multipole(3))
671 : ! Obtain dipoles in MO basis
672 : CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
673 : qs_env, mo_coeff, ref_point_multipole, 1, &
674 42 : homo, virtual, fm_eigvec_X%matrix_struct%context)
675 : ! Compute exciton descriptors from these multipoles
676 42 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
677 : ! Obtain quadrupoles in MO basis
678 40 : ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
679 40 : ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
680 40 : ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
681 : CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
682 : qs_env, mo_coeff, ref_point_multipole, 2, &
683 4 : homo, virtual, fm_eigvec_X%matrix_struct%context)
684 : ! Iterate over excitation index outside of routine to make it compatible with tddft module
685 236 : ALLOCATE (exc_descr(mp2_env%bse%num_print_exc_descr))
686 104 : DO i_exc = 1, mp2_env%bse%num_print_exc_descr
687 : CALL reshuffle_eigvec(fm_eigvec_X, fm_X_ia, homo, virtual, i_exc, &
688 100 : .FALSE., unit_nr, mp2_env)
689 100 : IF (.NOT. flag_TDA) THEN
690 : CALL reshuffle_eigvec(fm_eigvec_Y, fm_Y_ia, homo, virtual, i_exc, &
691 50 : .FALSE., unit_nr, mp2_env)
692 :
693 : CALL get_exciton_descriptors(exc_descr, fm_X_ia, &
694 : fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
695 : fm_quadpole_ai_trunc, &
696 : i_exc, homo, virtual, &
697 50 : fm_Y_ia)
698 : ELSE
699 : CALL get_exciton_descriptors(exc_descr, fm_X_ia, &
700 : fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
701 : fm_quadpole_ai_trunc, &
702 50 : i_exc, homo, virtual)
703 : END IF
704 100 : CALL cp_fm_release(fm_X_ia)
705 104 : IF (.NOT. flag_TDA) THEN
706 50 : CALL cp_fm_release(fm_Y_ia)
707 : END IF
708 : END DO
709 : END IF
710 :
711 42 : IF (mp2_env%bse%bse_spin_config == 0) THEN
712 : CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
713 : trans_mom_bse, oscill_str, polarizability_residues, &
714 : mp2_env, homo, virtual, unit_nr, &
715 42 : fm_eigvec_Y)
716 : END IF
717 :
718 : ! Prints basic definitions used in BSE calculation
719 : CALL print_output_header(homo, virtual, homo_irred, flag_TDA, &
720 42 : multiplet, alpha, mp2_env, unit_nr)
721 :
722 : ! Prints excitation energies up to user-specified number
723 : CALL print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
724 42 : info_approximation, mp2_env, unit_nr)
725 :
726 : ! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
727 : CALL print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
728 42 : info_approximation, mp2_env, unit_nr, fm_eigvec_Y)
729 :
730 : ! Prints optical properties, if state is a singlet
731 : CALL print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
732 : homo, virtual, homo_irred, flag_TDA, &
733 42 : info_approximation, mp2_env, unit_nr)
734 : ! Print exciton descriptors if keyword is invoked
735 42 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
736 : CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
737 : mp2_env%bse%num_print_exc_descr, mp2_env%bse%bse_debug_print, &
738 : mp2_env%bse%print_directional_exc_descr, &
739 4 : 'BSE|', qs_env)
740 : END IF
741 :
742 : ! Compute and print excitation wavefunctions
743 42 : IF (mp2_env%bse%do_nto_analysis) THEN
744 4 : IF (unit_nr > 0) THEN
745 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
746 : WRITE (unit_nr, '(T2,A4,T7,A47)') &
747 2 : 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
748 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
749 : END IF
750 : CALL calculate_NTOs(fm_eigvec_X, fm_eigvec_Y, &
751 : mo_coeff, homo, virtual, &
752 : info_approximation, &
753 : oscill_str, &
754 4 : qs_env, unit_nr, mp2_env)
755 : END IF
756 :
757 168 : DO idir = 1, n_moments_di
758 126 : CALL cp_fm_release(fm_dipole_ai_trunc(idir))
759 126 : CALL cp_fm_release(fm_dipole_ij_trunc(idir))
760 168 : CALL cp_fm_release(fm_dipole_ab_trunc(idir))
761 : END DO
762 42 : IF (mp2_env%bse%num_print_exc_descr > 0) THEN
763 40 : DO idir = 1, n_moments_quad
764 36 : CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
765 36 : CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
766 40 : CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
767 : END DO
768 4 : DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
769 4 : DEALLOCATE (exc_descr)
770 : END IF
771 42 : DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
772 42 : DEALLOCATE (ref_point_multipole)
773 42 : IF (mp2_env%bse%bse_spin_config == 0) THEN
774 42 : DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
775 : END IF
776 42 : IF (unit_nr > 0) THEN
777 21 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
778 21 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
779 : END IF
780 :
781 42 : CALL timestop(handle)
782 :
783 84 : END SUBROUTINE postprocess_bse
784 :
785 : END MODULE bse_full_diag
|
