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_util, ONLY: comp_eigvec_coeff_BSE,&
17 : compute_absorption_spectrum,&
18 : filter_eigvec_contrib,&
19 : fm_general_add_bse,&
20 : get_oscillator_strengths
21 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
22 : cp_blacs_env_release,&
23 : cp_blacs_env_type
24 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
25 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
26 : cp_fm_power
27 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_get_info,&
32 : cp_fm_release,&
33 : cp_fm_set_all,&
34 : cp_fm_to_fm,&
35 : cp_fm_type
36 : USE input_constants, ONLY: bse_singlet,&
37 : bse_triplet
38 : USE kinds, ONLY: dp
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE mp2_types, ONLY: mp2_type
41 : USE parallel_gemm_api, ONLY: parallel_gemm
42 : USE physcon, ONLY: evolt
43 : USE qs_environment_types, ONLY: qs_environment_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : PRIVATE
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag'
51 :
52 : PUBLIC :: create_A, diagonalize_A, create_B, create_hermitian_form_of_ABBA, &
53 : diagonalize_C
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
59 : !> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
60 : !> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
61 : !> α is a spin-dependent factor
62 : !> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
63 : !> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
64 : !> \param fm_mat_S_ia_bse ...
65 : !> \param fm_mat_S_bar_ij_bse ...
66 : !> \param fm_mat_S_ab_bse ...
67 : !> \param fm_A ...
68 : !> \param Eigenval ...
69 : !> \param unit_nr ...
70 : !> \param homo ...
71 : !> \param virtual ...
72 : !> \param dimen_RI ...
73 : !> \param mp2_env ...
74 : !> \param para_env ...
75 : ! **************************************************************************************************
76 18 : SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
77 18 : fm_A, Eigenval, unit_nr, &
78 : homo, virtual, dimen_RI, mp2_env, &
79 : para_env)
80 :
81 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, &
82 : fm_mat_S_ab_bse
83 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
84 : REAL(KIND=dp), DIMENSION(:) :: Eigenval
85 : INTEGER, INTENT(IN) :: unit_nr, homo, virtual, dimen_RI
86 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
87 : TYPE(mp_para_env_type), INTENT(INOUT) :: para_env
88 :
89 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_A'
90 :
91 : INTEGER :: a_virt_row, handle, i_occ_row, &
92 : i_row_global, ii, j_col_global, jj, &
93 : ncol_local_A, nrow_local_A
94 : INTEGER, DIMENSION(4) :: reordering
95 18 : INTEGER, DIMENSION(:), POINTER :: col_indices_A, row_indices_A
96 : REAL(KIND=dp) :: alpha, eigen_diff
97 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
98 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_A, fm_struct_W
99 : TYPE(cp_fm_type) :: fm_W
100 :
101 18 : CALL timeset(routineN, handle)
102 :
103 18 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
104 0 : WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
105 : END IF
106 :
107 : !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
108 36 : SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
109 : CASE (bse_singlet)
110 18 : alpha = 2.0_dp
111 : CASE (bse_triplet)
112 18 : alpha = 0.0_dp
113 : END SELECT
114 :
115 : ! create the blacs env for ij matrices (NOT fm_mat_S_ia_bse%matrix_struct related parallel_gemms!)
116 18 : NULLIFY (blacs_env)
117 18 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
118 :
119 : ! We have to use the same blacs_env for A as for the matrices fm_mat_S_ia_bse from RPA
120 : ! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
121 : ! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab
122 : ! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb
123 : ! We use the A matrix already from the start instead of v
124 : CALL cp_fm_struct_create(fm_struct_A, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
125 18 : ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
126 18 : CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb")
127 18 : CALL cp_fm_set_all(fm_A, 0.0_dp)
128 :
129 : CALL cp_fm_struct_create(fm_struct_W, context=fm_mat_S_ab_bse%matrix_struct%context, nrow_global=homo**2, &
130 18 : ncol_global=virtual**2, para_env=fm_mat_S_ab_bse%matrix_struct%para_env)
131 18 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab")
132 18 : CALL cp_fm_set_all(fm_W, 0.0_dp)
133 :
134 : ! Create A matrix from GW Energies, v_ia,jb and W_ij,ab (different blacs_env!)
135 : ! v_ia,jb, which is directly initialized in A (with a factor of alpha)
136 : ! v_ia,jb = \sum_P B^P_ia B^P_jb
137 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
138 : matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
139 18 : matrix_c=fm_A)
140 :
141 18 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
142 0 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
143 : END IF
144 :
145 : !W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab
146 : CALL parallel_gemm(transa="T", transb="N", m=homo**2, n=virtual**2, k=dimen_RI, alpha=1.0_dp, &
147 : matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, &
148 18 : matrix_c=fm_W)
149 :
150 18 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
151 0 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
152 : END IF
153 :
154 : ! We start by moving data from local parts of W_ij,ab to the full matrix A_ia,jb using buffers
155 : CALL cp_fm_get_info(matrix=fm_A, &
156 : nrow_local=nrow_local_A, &
157 : ncol_local=ncol_local_A, &
158 : row_indices=row_indices_A, &
159 18 : col_indices=col_indices_A)
160 : ! Writing -1.0_dp * W_ij,ab to A_ia,jb, i.e. beta = -1.0_dp,
161 : ! W_ij,ab: nrow_secidx_in = homo, ncol_secidx_in = virtual
162 : ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
163 18 : reordering = (/1, 3, 2, 4/)
164 : CALL fm_general_add_bse(fm_A, fm_W, -1.0_dp, homo, virtual, &
165 18 : virtual, virtual, unit_nr, reordering, mp2_env)
166 : !full matrix W is not needed anymore, release it to save memory
167 18 : CALL cp_fm_release(fm_W)
168 :
169 : !Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb
170 450 : DO ii = 1, nrow_local_A
171 :
172 432 : i_row_global = row_indices_A(ii)
173 :
174 21186 : DO jj = 1, ncol_local_A
175 :
176 20736 : j_col_global = col_indices_A(jj)
177 :
178 21168 : IF (i_row_global == j_col_global) THEN
179 432 : i_occ_row = (i_row_global - 1)/virtual + 1
180 432 : a_virt_row = MOD(i_row_global - 1, virtual) + 1
181 432 : eigen_diff = Eigenval(a_virt_row + homo) - Eigenval(i_occ_row)
182 432 : fm_A%local_data(ii, jj) = fm_A%local_data(ii, jj) + eigen_diff
183 :
184 : END IF
185 : END DO
186 : END DO
187 :
188 18 : CALL cp_fm_struct_release(fm_struct_A)
189 18 : CALL cp_fm_struct_release(fm_struct_W)
190 :
191 18 : CALL cp_blacs_env_release(blacs_env)
192 :
193 18 : CALL timestop(handle)
194 :
195 72 : END SUBROUTINE create_A
196 :
197 : ! **************************************************************************************************
198 : !> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
199 : !> B_ia,jb = α * v_ia,jb - W_ib,aj
200 : !> α is a spin-dependent factor
201 : !> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
202 : !> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
203 : !> \param fm_mat_S_ia_bse ...
204 : !> \param fm_mat_S_bar_ia_bse ...
205 : !> \param fm_B ...
206 : !> \param homo ...
207 : !> \param virtual ...
208 : !> \param dimen_RI ...
209 : !> \param unit_nr ...
210 : !> \param mp2_env ...
211 : ! **************************************************************************************************
212 6 : SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
213 : homo, virtual, dimen_RI, unit_nr, mp2_env)
214 :
215 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse
216 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_B
217 : INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr
218 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
219 :
220 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_B'
221 :
222 : INTEGER :: handle
223 : INTEGER, DIMENSION(4) :: reordering
224 : REAL(KIND=dp) :: alpha
225 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_v
226 : TYPE(cp_fm_type) :: fm_W
227 :
228 6 : CALL timeset(routineN, handle)
229 :
230 6 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
231 0 : WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
232 : END IF
233 :
234 : ! Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
235 12 : SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
236 : CASE (bse_singlet)
237 6 : alpha = 2.0_dp
238 : CASE (bse_triplet)
239 6 : alpha = 0.0_dp
240 : END SELECT
241 :
242 : CALL cp_fm_struct_create(fm_struct_v, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
243 6 : ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
244 6 : CALL cp_fm_create(fm_B, fm_struct_v, name="fm_B_iajb")
245 6 : CALL cp_fm_set_all(fm_B, 0.0_dp)
246 :
247 6 : CALL cp_fm_create(fm_W, fm_struct_v, name="fm_W_ibaj")
248 6 : CALL cp_fm_set_all(fm_W, 0.0_dp)
249 :
250 6 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
251 0 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
252 : END IF
253 : ! v_ia,jb = \sum_P B^P_ia B^P_jb
254 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
255 : matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
256 6 : matrix_c=fm_B)
257 :
258 : ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj
259 : CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=1.0_dp, &
260 : matrix_a=fm_mat_S_bar_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
261 6 : matrix_c=fm_W)
262 : ! from W_ib,ja to A_ia,jb (formally: W_ib,aj, but our internal indexorder is different)
263 : ! Writing -1.0_dp * W_ib,ja to A_ia,jb, i.e. beta = -1.0_dp,
264 : ! W_ib,ja: nrow_secidx_in = virtual, ncol_secidx_in = virtual
265 : ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
266 6 : reordering = (/1, 4, 3, 2/)
267 : CALL fm_general_add_bse(fm_B, fm_W, -1.0_dp, virtual, virtual, &
268 6 : virtual, virtual, unit_nr, reordering, mp2_env)
269 :
270 6 : CALL cp_fm_release(fm_W)
271 6 : CALL cp_fm_struct_release(fm_struct_v)
272 6 : CALL timestop(handle)
273 :
274 18 : END SUBROUTINE create_B
275 :
276 : ! **************************************************************************************************
277 : !> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
278 : !> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
279 : !> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
280 : !> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
281 : !> \param fm_A ...
282 : !> \param fm_B ...
283 : !> \param fm_C ...
284 : !> \param fm_sqrt_A_minus_B ...
285 : !> \param fm_inv_sqrt_A_minus_B ...
286 : !> \param homo ...
287 : !> \param virtual ...
288 : !> \param unit_nr ...
289 : !> \param mp2_env ...
290 : !> \param diag_est ...
291 : ! **************************************************************************************************
292 48 : SUBROUTINE create_hermitian_form_of_ABBA(fm_A, fm_B, fm_C, &
293 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
294 : homo, virtual, unit_nr, mp2_env, diag_est)
295 :
296 : TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B
297 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C, fm_sqrt_A_minus_B, &
298 : fm_inv_sqrt_A_minus_B
299 : INTEGER, INTENT(IN) :: homo, virtual, unit_nr
300 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
301 : REAL(KIND=dp), INTENT(IN) :: diag_est
302 :
303 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_hermitian_form_of_ABBA'
304 :
305 : INTEGER :: dim_mat, handle, n_dependent
306 : REAL(KIND=dp), DIMENSION(2) :: eigvals_AB_diff
307 : TYPE(cp_fm_type) :: fm_A_minus_B, fm_A_plus_B, fm_dummy, &
308 : fm_work_product
309 :
310 6 : CALL timeset(routineN, handle)
311 :
312 6 : IF (unit_nr > 0) THEN
313 3 : WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
314 6 : ' with size of A. This will take around ', diag_est, " s."
315 : END IF
316 :
317 : ! Create work matrices, which will hold A+B and A-B and their powers
318 : ! C is created afterwards to save memory
319 : ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
320 : ! \_______/ \___/ \______/
321 : ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
322 : ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
323 : ! Intermediate work matrices:
324 : ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
325 : ! fm_A_minus_B: (A-B) EQ.III
326 : ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
327 6 : CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct)
328 6 : CALL cp_fm_to_fm(fm_A, fm_A_plus_B)
329 6 : CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct)
330 6 : CALL cp_fm_to_fm(fm_A, fm_A_minus_B)
331 6 : CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct)
332 6 : CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp)
333 6 : CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct)
334 6 : CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp)
335 :
336 6 : CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
337 :
338 6 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
339 0 : WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
340 : END IF
341 :
342 : ! Add/Substract B (cf. EQs. Ib and III)
343 6 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B)
344 6 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_minus_B, -1.0_dp, fm_B)
345 :
346 : ! cp_fm_power will overwrite matrix, therefore we create copies
347 6 : CALL cp_fm_to_fm(fm_A_minus_B, fm_inv_sqrt_A_minus_B)
348 :
349 : ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
350 : ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
351 :
352 : ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
353 6 : CALL cp_fm_create(fm_dummy, fm_A%matrix_struct)
354 : ! Create (A-B)^-0.5 (cf. EQ.II)
355 6 : CALL cp_fm_power(fm_inv_sqrt_A_minus_B, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_AB_diff)
356 6 : CALL cp_fm_release(fm_dummy)
357 : ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
358 : ! In this case, the procedure for hermitian form of ABBA is not applicable
359 6 : IF (eigvals_AB_diff(1) < 0) THEN
360 : CALL cp_abort(__LOCATION__, &
361 : "Matrix (A-B) is not positive definite. "// &
362 0 : "Hermitian diagonalization of full ABBA matrix is ill-defined.")
363 : END IF
364 :
365 : ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
366 : ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
367 : ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
368 6 : dim_mat = homo*virtual
369 : 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, &
370 6 : fm_sqrt_A_minus_B)
371 :
372 : ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
373 : 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, &
374 6 : fm_work_product)
375 :
376 : ! Release to save memory
377 6 : CALL cp_fm_release(fm_A_plus_B)
378 6 : CALL cp_fm_release(fm_A_minus_B)
379 :
380 : ! Now create full
381 6 : CALL cp_fm_create(fm_C, fm_A%matrix_struct)
382 6 : CALL cp_fm_set_all(fm_C, 0.0_dp)
383 : ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
384 : 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, &
385 6 : fm_C)
386 6 : CALL cp_fm_release(fm_work_product)
387 :
388 6 : IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
389 0 : WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
390 : END IF
391 :
392 6 : CALL timestop(handle)
393 6 : END SUBROUTINE
394 :
395 : ! **************************************************************************************************
396 : !> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
397 : !> Here, the eigenvectors Z^n relate to X^n via
398 : !> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
399 : !> \param fm_C ...
400 : !> \param homo ...
401 : !> \param virtual ...
402 : !> \param homo_irred ...
403 : !> \param fm_sqrt_A_minus_B ...
404 : !> \param fm_inv_sqrt_A_minus_B ...
405 : !> \param unit_nr ...
406 : !> \param diag_est ...
407 : !> \param mp2_env ...
408 : !> \param qs_env ...
409 : !> \param mo_coeff ...
410 : ! **************************************************************************************************
411 6 : SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
412 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
413 6 : unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
414 :
415 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C
416 : INTEGER, INTENT(IN) :: homo, virtual, homo_irred
417 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B
418 : INTEGER, INTENT(IN) :: unit_nr
419 : REAL(KIND=dp), INTENT(IN) :: diag_est
420 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
421 : TYPE(qs_environment_type), POINTER :: qs_env
422 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
423 :
424 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_C'
425 :
426 : INTEGER :: diag_info, handle
427 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
428 : TYPE(cp_fm_type) :: fm_eigvec_X, fm_eigvec_Y, fm_eigvec_Z, &
429 : fm_mat_eigvec_transform_diff, &
430 : fm_mat_eigvec_transform_sum
431 :
432 6 : CALL timeset(routineN, handle)
433 :
434 6 : IF (unit_nr > 0) THEN
435 3 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
436 6 : 'This will take around ', diag_est, ' s.'
437 : END IF
438 :
439 : !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
440 : !Now: Diagonalize it
441 6 : CALL cp_fm_create(fm_eigvec_Z, fm_C%matrix_struct)
442 :
443 18 : ALLOCATE (Exc_ens(homo*virtual))
444 :
445 6 : CALL choose_eigv_solver(fm_C, fm_eigvec_Z, Exc_ens, diag_info)
446 6 : CPASSERT(diag_info == 0)
447 : ! C could have negative eigenvalues, since we do not explicitly check A+B
448 : ! for positive definiteness (would make another O(N^6) Diagon. necessary)
449 : ! Instead, we include a check here
450 6 : IF (Exc_ens(1) < 0) THEN
451 : CALL cp_abort(__LOCATION__, &
452 : "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
453 0 : "(A+B) is not positive definite.")
454 : END IF
455 294 : Exc_ens = SQRT(Exc_ens)
456 :
457 : ! Prepare eigenvector for interpretation of singleparticle transitions
458 : ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
459 : ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
460 :
461 : ! Following Furche, we basically use Eqs. (A10): First, we multiply
462 : ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
463 : ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
464 :
465 : ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
466 6 : CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_C%matrix_struct)
467 6 : CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
468 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
469 : matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
470 6 : matrix_c=fm_mat_eigvec_transform_sum)
471 6 : CALL cp_fm_release(fm_sqrt_A_minus_B)
472 6 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_sum, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
473 :
474 : ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
475 6 : CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_C%matrix_struct)
476 6 : CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
477 : CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
478 : matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
479 6 : matrix_c=fm_mat_eigvec_transform_diff)
480 6 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
481 6 : CALL cp_fm_release(fm_eigvec_Z)
482 6 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_diff, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
483 :
484 : ! Now, we add the two equations to obtain X_n
485 : ! Add overwrites the first argument, therefore we copy it beforehand
486 6 : CALL cp_fm_create(fm_eigvec_X, fm_C%matrix_struct)
487 6 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_X)
488 6 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_X, 1.0_dp, fm_mat_eigvec_transform_diff)
489 :
490 : ! Now, we subtract the two equations to obtain Y_n
491 : ! Add overwrites the first argument, therefore we copy it beforehand
492 6 : CALL cp_fm_create(fm_eigvec_Y, fm_C%matrix_struct)
493 6 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_Y)
494 6 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_Y, -1.0_dp, fm_mat_eigvec_transform_diff)
495 :
496 : !Cleanup
497 6 : CALL cp_fm_release(fm_mat_eigvec_transform_diff)
498 6 : CALL cp_fm_release(fm_mat_eigvec_transform_sum)
499 :
500 : CALL success_message(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
501 : homo, virtual, homo_irred, unit_nr, &
502 6 : .FALSE., fm_eigvec_Y)
503 :
504 6 : DEALLOCATE (Exc_ens)
505 6 : CALL cp_fm_release(fm_eigvec_X)
506 6 : CALL cp_fm_release(fm_eigvec_Y)
507 :
508 6 : CALL timestop(handle)
509 :
510 42 : END SUBROUTINE diagonalize_C
511 :
512 : ! **************************************************************************************************
513 : !> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
514 : !> \param fm_A ...
515 : !> \param homo ...
516 : !> \param virtual ...
517 : !> \param homo_irred ...
518 : !> \param unit_nr ...
519 : !> \param diag_est ...
520 : !> \param mp2_env ...
521 : !> \param qs_env ...
522 : !> \param mo_coeff ...
523 : ! **************************************************************************************************
524 12 : SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
525 12 : unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
526 :
527 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
528 : INTEGER, INTENT(IN) :: homo, virtual, homo_irred, unit_nr
529 : REAL(KIND=dp), INTENT(IN) :: diag_est
530 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
531 : TYPE(qs_environment_type), POINTER :: qs_env
532 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
533 :
534 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A'
535 :
536 : INTEGER :: diag_info, handle
537 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
538 : TYPE(cp_fm_type) :: fm_eigvec
539 :
540 12 : CALL timeset(routineN, handle)
541 :
542 : !Continue with formatting of subroutine create_A
543 12 : IF (unit_nr > 0) THEN
544 6 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
545 12 : 'This will take around ', diag_est, ' s.'
546 : END IF
547 :
548 : !We have now the full matrix A_iajb, distributed over all ranks
549 : !Now: Diagonalize it
550 12 : CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
551 :
552 36 : ALLOCATE (Exc_ens(homo*virtual))
553 :
554 12 : CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
555 12 : CPASSERT(diag_info == 0)
556 : CALL success_message(Exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
557 12 : homo, virtual, homo_irred, unit_nr, .TRUE.)
558 :
559 12 : CALL cp_fm_release(fm_eigvec)
560 12 : DEALLOCATE (Exc_ens)
561 :
562 12 : CALL timestop(handle)
563 :
564 36 : END SUBROUTINE diagonalize_A
565 :
566 : ! **************************************************************************************************
567 : !> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
568 : !> \param Exc_ens ...
569 : !> \param fm_eigvec_X ...
570 : !> \param mp2_env ...
571 : !> \param qs_env ...
572 : !> \param mo_coeff ...
573 : !> \param homo ...
574 : !> \param virtual ...
575 : !> \param homo_irred ...
576 : !> \param unit_nr ...
577 : !> \param flag_TDA ...
578 : !> \param fm_eigvec_Y ...
579 : ! **************************************************************************************************
580 18 : SUBROUTINE success_message(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
581 : homo, virtual, homo_irred, unit_nr, &
582 : flag_TDA, fm_eigvec_Y)
583 :
584 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
585 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
586 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
587 : TYPE(qs_environment_type), POINTER :: qs_env
588 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
589 : INTEGER :: homo, virtual, homo_irred, unit_nr
590 : LOGICAL, OPTIONAL :: flag_TDA
591 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
592 :
593 : CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message'
594 :
595 : CHARACTER(LEN=10) :: info_approximation, multiplet
596 : INTEGER :: handle, i_exc, k, num_entries
597 18 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
598 : REAL(KIND=dp) :: alpha
599 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries, oscill_str
600 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dipole_exc
601 :
602 18 : CALL timeset(routineN, handle)
603 :
604 : !Prepare variables for printing
605 18 : IF (mp2_env%ri_g0w0%bse_spin_config == 0) THEN
606 18 : multiplet = "Singlet"
607 18 : alpha = 2.0_dp
608 : ELSE
609 0 : multiplet = "Triplet"
610 0 : alpha = 0.0_dp
611 : END IF
612 18 : IF (.NOT. PRESENT(flag_TDA)) THEN
613 0 : flag_TDA = .FALSE.
614 : END IF
615 18 : IF (flag_TDA) THEN
616 12 : info_approximation = " -TDA- "
617 : ELSE
618 6 : info_approximation = "-ABBA-"
619 : END IF
620 :
621 : !Get information about eigenvector
622 :
623 18 : IF (unit_nr > 0) THEN
624 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
625 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
626 9 : IF (flag_TDA) THEN
627 6 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
628 6 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *'
629 6 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
630 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
631 6 : WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
632 12 : 'the BSE within the TDA:'
633 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
634 6 : WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
635 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
636 6 : WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
637 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
638 6 : WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb X_jb^n ) = Ω^n X_ia^n'
639 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
640 6 : WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
641 6 : WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
642 : ELSE
643 3 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
644 3 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full ("ABBA") Bethe Salpeter equation (BSE) (i.e. without TDA) *'
645 3 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
646 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
647 3 : WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
648 6 : 'the BSE without the TDA:'
649 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
650 3 : WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|'
651 3 : WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
652 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
653 3 : WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
654 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
655 3 : WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', ' sum_jb ( A_ia,jb X_jb^n + B_ia,jb Y_jb^n ) = Ω^n X_ia^n'
656 3 : WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', '- sum_jb ( B_ia,jb X_jb^n + A_ia,jb Y_jb^n ) = Ω^n Y_ia^n'
657 : END IF
658 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
659 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
660 9 : WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
661 18 : 'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
662 9 : WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
663 18 : 'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
664 9 : WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
665 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
666 9 : WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
667 9 : IF (.NOT. flag_TDA) THEN
668 3 : WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
669 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
670 3 : WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
671 3 : WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
672 : END IF
673 9 : IF (.NOT. flag_TDA) THEN
674 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
675 3 : WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', 'The BSE is solved for Ω^n and X_ia^n as a hermitian problem, e.g. Eq.(42)'
676 3 : WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
677 : ! WRITE (unit_nr, '(T2,A4,T7,A69)') 'BSE|', 'C_ia,jb = sum_kc,ld ((A-B)^0.5)_ia,kc (A+B)_kc,ld ((A-B)^0.5)_ld,jb'
678 : ! WRITE (unit_nr, '(T2,A4)') 'BSE|'
679 : ! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X+Y)_ia,n = sum_jb,m ((A-B)^0.5)ia,jb Z_jb,m (Ω_m)^-0.5'
680 : ! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X-Y)_ia,n = sum_jb,m ((A-B)^-0.5)ia,jb Z_jb,m (Ω_m)^0.5'
681 : END IF
682 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
683 9 : WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
684 9 : WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
685 9 : WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
686 9 : WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
687 9 : WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction'
688 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
689 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
690 9 : WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
691 18 : 'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
692 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
693 9 : IF (flag_TDA) THEN
694 6 : WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
695 : ELSE
696 3 : WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
697 : END IF
698 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
699 9 : WRITE (unit_nr, '(T2,A4,T11,A12,T26,A11,T44,A8,T55,A27)') 'BSE|', &
700 18 : 'Excitation n', "Spin Config", 'TDA/ABBA', 'Excitation energy Ω^n (eV)'
701 : END IF
702 : !prints actual energies values
703 18 : IF (unit_nr > 0) THEN
704 234 : DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
705 : WRITE (unit_nr, '(T2,A4,T7,I16,T30,A7,T46,A6,T59,F22.4)') &
706 234 : 'BSE|', i_exc, multiplet, info_approximation, Exc_ens(i_exc)*evolt
707 : END DO
708 : END IF
709 : !prints single particle transitions
710 18 : IF (unit_nr > 0) THEN
711 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
712 :
713 : WRITE (unit_nr, '(T2,A4,T7,A70)') &
714 9 : 'BSE|', "Excitations are built up by the following single-particle transitions,"
715 : WRITE (unit_nr, '(T2,A4,T7,A42,F5.2,A2)') &
716 9 : 'BSE|', "neglecting contributions where |X_ia^n| < ", mp2_env%ri_g0w0%eps_x, " :"
717 :
718 9 : WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
719 18 : homo_irred, ' and LUMO a =', homo_irred + 1, " --"
720 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
721 : WRITE (unit_nr, '(T2,A4,T7,A12,T25,A1,T27,A2,T34,A1,T54,A8,T71,A10)') &
722 9 : "BSE|", "Excitation n", "i", "=>", "a", 'TDA/ABBA', "|X_ia^n|"
723 : END IF
724 468 : DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
725 : !Iterate through eigenvector and print values above threshold
726 : CALL filter_eigvec_contrib(fm_eigvec_X, idx_homo, idx_virt, eigvec_entries, &
727 450 : i_exc, virtual, num_entries, mp2_env)
728 450 : IF (unit_nr > 0) THEN
729 225 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
730 756 : DO k = 1, num_entries
731 : WRITE (unit_nr, '(T2,A4,T14,I5,T21,I5,T27,A2,T30,I5,T56,A6,T65,F16.4)') &
732 531 : "BSE|", i_exc, homo_irred - homo + idx_homo(k), "=>", &
733 1287 : homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k))
734 : END DO
735 : END IF
736 :
737 450 : DEALLOCATE (idx_homo)
738 450 : DEALLOCATE (idx_virt)
739 468 : DEALLOCATE (eigvec_entries)
740 : END DO
741 : ! Compute BSE dipoles and oscillator strengths
742 : CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, &
743 : dipole_exc, oscill_str, &
744 : mp2_env, qs_env, mo_coeff, &
745 : homo, virtual, unit_nr, &
746 18 : fm_eigvec_Y)
747 18 : IF (unit_nr > 0) THEN
748 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
749 : WRITE (unit_nr, '(T2,A4,T7,A63)') &
750 9 : 'BSE|', "Dipoles d_r^n and oscillator strength f^n of excitation level n"
751 : WRITE (unit_nr, '(T2,A4,T7,A17)') &
752 9 : 'BSE|', "are obtained from"
753 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
754 9 : IF (flag_TDA) THEN
755 : WRITE (unit_nr, '(T2,A4,T12,A51)') &
756 6 : 'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | µ | ψ_a > X_ia^n"
757 : ELSE
758 : WRITE (unit_nr, '(T2,A4,T12,A63)') &
759 3 : 'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | µ | ψ_a > ( X_ia^n + Y_ia^n )"
760 : END IF
761 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
762 : WRITE (unit_nr, '(T2,A4,T20,A44)') &
763 9 : 'BSE|', "f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2"
764 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
765 : WRITE (unit_nr, '(T2,A4,T7,A19)') &
766 9 : 'BSE|', "where we introduced"
767 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
768 : WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
769 9 : 'BSE|', "ψ_i:", "occupied molecular orbitals,"
770 : WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
771 9 : 'BSE|', "ψ_a:", "empty molecular orbitals and"
772 : WRITE (unit_nr, '(T2,A4,T7,A3,T15,A16)') &
773 9 : 'BSE|', "µ:", "dipole operator."
774 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
775 : WRITE (unit_nr, '(T2,A4,T7,A28)') &
776 9 : 'BSE|', "prelim Ref.: Eqs. (23), (24)"
777 : WRITE (unit_nr, '(T2,A4,T7,A71)') &
778 9 : 'BSE|', "in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290"
779 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
780 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
781 9 : WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T35,A8,T45,A8,T55,A8,T65,A16)') 'BSE|', &
782 18 : 'Excitation n', "TDA/ABBA", "Dipole X", "Dipole Y", "Dipole Z", 'Osc. strength'
783 234 : DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
784 : WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T35,F8.3,T45,F8.3,T55,F8.3,T65,F16.3)') &
785 225 : 'BSE|', i_exc, info_approximation, dipole_exc(1, 1, i_exc), dipole_exc(2, 1, i_exc), &
786 459 : dipole_exc(3, 1, i_exc), oscill_str(i_exc)
787 : END DO
788 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
789 : END IF
790 :
791 : ! Compute and print the absorption spectrum to external file
792 18 : IF (mp2_env%ri_g0w0%bse_print_spectrum) THEN
793 : CALL compute_absorption_spectrum(oscill_str, Exc_ens, &
794 0 : info_approximation, unit_nr, mp2_env)
795 : END IF
796 18 : DEALLOCATE (oscill_str, dipole_exc)
797 18 : IF (unit_nr > 0) THEN
798 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
799 9 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
800 : END IF
801 :
802 18 : CALL timestop(handle)
803 :
804 18 : END SUBROUTINE success_message
805 :
806 : END MODULE bse_full_diag
|
