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