(git:374b731)
Loading...
Searching...
No Matches
bse_full_diag.F
Go to the documentation of this file.
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! **************************************************************************************************
15
28 USE cp_fm_types, ONLY: cp_fm_create,&
34 USE input_constants, ONLY: bse_singlet,&
36 USE kinds, ONLY: dp
38 USE mp2_types, ONLY: mp2_type
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
51
52CONTAINS
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 SUBROUTINE create_a(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
74 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 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 CALL timeset(routinen, handle)
99
100 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
101 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 SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
106 CASE (bse_singlet)
107 alpha = 2.0_dp
108 CASE (bse_triplet)
109 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 NULLIFY (blacs_env)
114 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 ncol_global=homo*virtual, para_env=fm_mat_s_ia_bse%matrix_struct%para_env)
123 CALL cp_fm_create(fm_a, fm_struct_a, name="fm_A_iajb")
124 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 ncol_global=virtual**2, para_env=fm_mat_s_ab_bse%matrix_struct%para_env)
128 CALL cp_fm_create(fm_w, fm_struct_w, name="fm_W_ijab")
129 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 matrix_c=fm_a)
137
138 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
139 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 matrix_c=fm_w)
146
147 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
148 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 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 reordering = (/1, 3, 2, 4/)
161 CALL fm_general_add_bse(fm_a, fm_w, -1.0_dp, homo, virtual, &
162 virtual, virtual, unit_nr, reordering, mp2_env)
163 !full matrix W is not needed anymore, release it to save memory
164 CALL cp_fm_struct_release(fm_struct_w)
165 CALL cp_fm_release(fm_w)
166 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 DO ii = 1, nrow_local_a
170
171 i_row_global = row_indices_a(ii)
172
173 DO jj = 1, ncol_local_a
174
175 j_col_global = col_indices_a(jj)
176
177 IF (i_row_global == j_col_global) THEN
178 i_occ_row = (i_row_global - 1)/virtual + 1
179 a_virt_row = mod(i_row_global - 1, virtual) + 1
180 eigen_diff = eigenval(a_virt_row + homo) - eigenval(i_occ_row)
181 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 CALL cp_fm_struct_release(fm_struct_a)
188 CALL cp_fm_struct_release(fm_struct_w)
189
190 CALL cp_blacs_env_release(blacs_env)
191
192 CALL timestop(handle)
193
194 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 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 CALL timeset(routinen, handle)
228
229 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
230 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 SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
235 CASE (bse_singlet)
236 alpha = 2.0_dp
237 CASE (bse_triplet)
238 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 ncol_global=homo*virtual, para_env=fm_mat_s_ia_bse%matrix_struct%para_env)
243 CALL cp_fm_create(fm_b, fm_struct_v, name="fm_B_iajb")
244 CALL cp_fm_set_all(fm_b, 0.0_dp)
245
246 CALL cp_fm_create(fm_w, fm_struct_v, name="fm_W_ibaj")
247 CALL cp_fm_set_all(fm_w, 0.0_dp)
248
249 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
250 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 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 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 reordering = (/1, 4, 3, 2/)
266 CALL fm_general_add_bse(fm_b, fm_w, -1.0_dp, virtual, virtual, &
267 virtual, virtual, unit_nr, reordering, mp2_env)
268
269 CALL cp_fm_release(fm_w)
270 CALL cp_fm_struct_release(fm_struct_v)
271 CALL timestop(handle)
272
273 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 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 CALL timeset(routinen, handle)
310
311 IF (unit_nr > 0) THEN
312 WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
313 ' 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 CALL cp_fm_create(fm_a_plus_b, fm_a%matrix_struct)
327 CALL cp_fm_to_fm(fm_a, fm_a_plus_b)
328 CALL cp_fm_create(fm_a_minus_b, fm_a%matrix_struct)
329 CALL cp_fm_to_fm(fm_a, fm_a_minus_b)
330 CALL cp_fm_create(fm_sqrt_a_minus_b, fm_a%matrix_struct)
331 CALL cp_fm_set_all(fm_sqrt_a_minus_b, 0.0_dp)
332 CALL cp_fm_create(fm_inv_sqrt_a_minus_b, fm_a%matrix_struct)
333 CALL cp_fm_set_all(fm_inv_sqrt_a_minus_b, 0.0_dp)
334
335 CALL cp_fm_create(fm_work_product, fm_a%matrix_struct)
336
337 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
338 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 CALL cp_fm_scale_and_add(1.0_dp, fm_a_plus_b, 1.0_dp, fm_b)
343 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 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 CALL cp_fm_create(fm_dummy, fm_a%matrix_struct)
353 ! Create (A-B)^-0.5 (cf. EQ.II)
354 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 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 IF (eigvals_ab_diff(1) < 0) THEN
359 CALL cp_abort(__location__, &
360 "Matrix (A-B) is not positive definite. "// &
361 "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 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 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 fm_work_product)
374
375 ! Release to save memory
376 CALL cp_fm_release(fm_a_plus_b)
377 CALL cp_fm_release(fm_a_minus_b)
378
379 ! Now create full
380 CALL cp_fm_create(fm_c, fm_a%matrix_struct)
381 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 fm_c)
385 CALL cp_fm_release(fm_work_product)
386
387 IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
388 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 CALL timestop(handle)
392 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 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 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 CALL timeset(routinen, handle)
427
428 IF (unit_nr > 0) THEN
429 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
430 '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 CALL cp_fm_create(fm_eigvec, fm_c%matrix_struct)
436
437 ALLOCATE (exc_ens(homo*virtual))
438
439 CALL choose_eigv_solver(fm_c, fm_eigvec, exc_ens, diag_info)
440 cpassert(diag_info == 0)
441 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 CALL cp_fm_create(fm_mat_eigvec_transform, fm_c%matrix_struct)
453 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 matrix_c=fm_mat_eigvec_transform)
457 CALL cp_fm_release(fm_sqrt_a_minus_b)
458 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 CALL cp_fm_create(fm_mat_eigvec_transform_neg, fm_c%matrix_struct)
462 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 matrix_c=fm_mat_eigvec_transform_neg)
466 CALL cp_fm_release(fm_inv_sqrt_a_minus_b)
467 CALL cp_fm_release(fm_eigvec)
468 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 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 CALL cp_fm_release(fm_mat_eigvec_transform_neg)
475
476 CALL success_message(exc_ens, fm_mat_eigvec_transform, mp2_env, &
477 homo, virtual, homo_irred, unit_nr, .false.)
478
479 DEALLOCATE (exc_ens)
480 CALL cp_fm_release(fm_mat_eigvec_transform)
481
482 CALL timestop(handle)
483
484 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
508 TYPE(cp_fm_type) :: fm_eigvec
509
510 CALL timeset(routinen, handle)
511
512 !Continue with formatting of subroutine create_A
513 IF (unit_nr > 0) THEN
514 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
515 '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 CALL cp_fm_create(fm_eigvec, fm_a%matrix_struct)
521
522 ALLOCATE (exc_ens(homo*virtual))
523
524 CALL choose_eigv_solver(fm_a, fm_eigvec, exc_ens, diag_info)
525 cpassert(diag_info == 0)
526 CALL success_message(exc_ens, fm_eigvec, mp2_env, &
527 homo, virtual, homo_irred, unit_nr, .true.)
528
529 CALL cp_fm_release(fm_eigvec)
530 DEALLOCATE (exc_ens)
531
532 CALL timestop(handle)
533
534 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
561 REAL(kind=dp) :: alpha
562 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
563
564 CALL timeset(routinen, handle)
565
566 !Prepare variables for printing
567 IF (mp2_env%ri_g0w0%bse_spin_config == 0) THEN
568 multiplet = "Singlet"
569 alpha = 2.0_dp
570 ELSE
571 multiplet = "Triplet"
572 alpha = 0.0_dp
573 END IF
574 IF (.NOT. PRESENT(flag_tda)) THEN
575 flag_tda = .false.
576 END IF
577 IF (flag_tda) THEN
578 info_approximation = " -TDA- "
579 ELSE
580 info_approximation = "-full-"
581 END IF
582
583 !Get information about eigenvector
584
585 IF (unit_nr > 0) THEN
586 WRITE (unit_nr, '(T2,A4)') 'BSE|'
587 WRITE (unit_nr, '(T2,A4)') 'BSE|'
588 IF (flag_tda) THEN
589 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
590 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *'
591 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
592 WRITE (unit_nr, '(T2,A4)') 'BSE|'
593 WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
594 'the BSE within the TDA:'
595 WRITE (unit_nr, '(T2,A4)') 'BSE|'
596 WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', Ω'A X^n = ^n X^n'
597 WRITE (unit_nr, '(T2,A4)') 'BSE|'
598 WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
599 WRITE (unit_nr, '(T2,A4)') 'BSE|'
600 WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', Ω'sum_jb ( A_ia,jb X_jb^n ) = ^n X_ia^n'
601 WRITE (unit_nr, '(T2,A4)') 'BSE|'
602 WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
603 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 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
606 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full Bethe Salpeter equation (BSE) (i.e. without TDA) *'
607 WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
608 WRITE (unit_nr, '(T2,A4)') 'BSE|'
609 WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
610 'the BSE without the TDA:'
611 WRITE (unit_nr, '(T2,A4)') 'BSE|'
612 WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|'
613 WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', Ω'|B A| |Y^n| = ^n |0 -1| |Y^n|'
614 WRITE (unit_nr, '(T2,A4)') 'BSE|'
615 WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
616 WRITE (unit_nr, '(T2,A4)') 'BSE|'
617 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 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 WRITE (unit_nr, '(T2,A4)') 'BSE|'
621 WRITE (unit_nr, '(T2,A4)') 'BSE|'
622 WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
623 'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
624 WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
625 'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
626 WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
627 WRITE (unit_nr, '(T2,A4)') 'BSE|'
628 WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', εεδδα'A_ia,jb = (_a-_i) _ij _ab + * v_ia,jb - W_ij,ab'
629 IF (.NOT. flag_tda) THEN
630 WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', α'B_ia,jb = * v_ia,jb - W_ib,aj'
631 WRITE (unit_nr, '(T2,A4)') 'BSE|'
632 WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
633 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 IF (.NOT. flag_tda) THEN
636 WRITE (unit_nr, '(T2,A4)') 'BSE|'
637 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 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 WRITE (unit_nr, '(T2,A4)') 'BSE|'
645 WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', ε'_...:', 'GW quasiparticle energy'
646 WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', δ'_...:', 'Kronecker delta'
647 WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', α':', 'spin-dependent factor (Singlet/Triplet)'
648 WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
649 WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction'
650 WRITE (unit_nr, '(T2,A4)') 'BSE|'
651 WRITE (unit_nr, '(T2,A4)') 'BSE|'
652 WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
653 'The spin-dependent factor is for the requested ', multiplet, α" is = ", alpha
654 WRITE (unit_nr, '(T2,A4)') 'BSE|'
655 IF (flag_tda) THEN
656 WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
657 ELSE
658 WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
659 END IF
660 WRITE (unit_nr, '(T2,A4)') 'BSE|'
661 WRITE (unit_nr, '(T2,A4,T13,A10,T27,A13,T42,A12,T59,A22)') 'BSE|', &
662 'Excitation', "Multiplet", 'TDA/full BSE', 'Excitation energy (eV)'
663 END IF
664 !prints actual energies values
665 IF (unit_nr > 0) THEN
666 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 'BSE|', i_exc, multiplet, " State", info_approximation, exc_ens(i_exc)*evolt
669 END DO
670 END IF
671 !prints single particle transitions
672 IF (unit_nr > 0) THEN
673 WRITE (unit_nr, '(T2,A4)') 'BSE|'
674
675 WRITE (unit_nr, '(T2,A4,T7,A70)') &
676 'BSE|', "Excitations are built up by the following single-particle transitions,"
677 WRITE (unit_nr, '(T2,A4,T7,A42,F5.2,A2)') &
678 'BSE|', "neglecting contributions where |X_ia^n| < ", mp2_env%ri_g0w0%eps_x, " :"
679
680 WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
681 homo_irred, ' and LUMO a =', homo_irred + 1, " --"
682 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 "BSE|", "n-th exc.", "i", "=>", "a", 'TDA/full BSE', "|X_ia^n|"
685 END IF
686 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 i_exc, virtual, num_entries, mp2_env)
690 IF (unit_nr > 0) THEN
691 WRITE (unit_nr, '(T2,A4)') 'BSE|'
692 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 "BSE|", i_exc, homo_irred - homo + idx_homo(k), "=>", &
695 homo_irred + idx_virt(k), info_approximation, abs(eigvec_entries(k))
696 END DO
697 END IF
698
699 DEALLOCATE (idx_homo)
700 DEALLOCATE (idx_virt)
701 DEALLOCATE (eigvec_entries)
702 END DO
703 IF (unit_nr > 0) THEN
704 WRITE (unit_nr, '(T2,A4)') 'BSE|'
705 WRITE (unit_nr, '(T2,A4)') 'BSE|'
706 END IF
707
708 CALL timestop(handle)
709 END SUBROUTINE success_message
710
711END MODULE bse_full_diag
Routines for the full diagonalization of GW + Bethe-Salpeter for computing electronic excitations.
subroutine, public diagonalize_c(fm_c, homo, virtual, homo_irred, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, diag_est, mp2_env)
Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . Here, the eigenvectors Z^n relate to X^n via Eq....
subroutine, public diagonalize_a(fm_a, homo, virtual, homo_irred, unit_nr, diag_est, mp2_env)
Solving hermitian eigenvalue equation A X^n = Ω^n X^n.
subroutine, public create_a(fm_mat_s_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_ab_bse, fm_a, eigenval, unit_nr, homo, virtual, dimen_ri, mp2_env, para_env)
Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W) A_ia,...
subroutine, public create_b(fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse, fm_b, homo, virtual, dimen_ri, unit_nr, mp2_env)
Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) B_ia,jb = α * v_ia,...
subroutine, public create_hermitian_form_of_abba(fm_a, fm_b, fm_c, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, homo, virtual, unit_nr, mp2_env, diag_est)
Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem (cf....
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
subroutine, public filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, i_exc, virtual, num_entries, mp2_env)
Filters eigenvector entries above a given threshold to describe excitations in the singleparticle bas...
Definition bse_util.F:1065
subroutine, public comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
Routine for computing the coefficients of the eigenvectors of the BSE matrix from a multiplication wi...
Definition bse_util.F:865
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Definition bse_util.F:151
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:896
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:208
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public bse_singlet
integer, parameter, public bse_triplet
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
Types needed for MP2 calculations.
Definition mp2_types.F:14
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
stores all the informations relevant to an mpi environment