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