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