Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Main routines for GW + Bethe-Salpeter for computing electronic excitations
10 : !> \par History
11 : !> 04.2024 created [Maximilian Graml]
12 : ! **************************************************************************************************
13 :
14 : MODULE bse_main
15 :
16 : USE bse_full_diag, ONLY: create_A,&
17 : create_B,&
18 : create_hermitian_form_of_ABBA,&
19 : diagonalize_A,&
20 : diagonalize_C
21 : USE bse_iterative, ONLY: do_subspace_iterations,&
22 : fill_local_3c_arrays
23 : USE bse_print, ONLY: print_BSE_start_flag
24 : USE bse_util, ONLY: adapt_BSE_input_params,&
25 : deallocate_matrices_bse,&
26 : estimate_BSE_resources,&
27 : mult_B_with_W,&
28 : truncate_BSE_matrices
29 : USE cp_control_types, ONLY: dft_control_type,&
30 : tddfpt2_control_type
31 : USE cp_fm_types, ONLY: cp_fm_release,&
32 : cp_fm_type
33 : USE cp_log_handling, ONLY: cp_get_default_logger,&
34 : cp_logger_type
35 : USE cp_output_handling, ONLY: debug_print_level
36 : USE group_dist_types, ONLY: group_dist_d1_type
37 : USE input_constants, ONLY: bse_abba,&
38 : bse_both,&
39 : bse_fulldiag,&
40 : bse_iterdiag,&
41 : bse_screening_alpha,&
42 : bse_screening_tdhf,&
43 : bse_tda
44 : USE kinds, ONLY: dp
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE mp2_types, ONLY: mp2_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
56 :
57 : PUBLIC :: start_bse_calculation
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Main subroutine managing BSE calculations
63 : !> \param fm_mat_S_ia_bse ...
64 : !> \param fm_mat_S_ij_bse ...
65 : !> \param fm_mat_S_ab_bse ...
66 : !> \param fm_mat_Q_static_bse_gemm ...
67 : !> \param Eigenval ...
68 : !> \param Eigenval_scf ...
69 : !> \param homo ...
70 : !> \param virtual ...
71 : !> \param dimen_RI ...
72 : !> \param dimen_RI_red ...
73 : !> \param bse_lev_virt ...
74 : !> \param gd_array ...
75 : !> \param color_sub ...
76 : !> \param mp2_env ...
77 : !> \param qs_env ...
78 : !> \param mo_coeff ...
79 : !> \param unit_nr ...
80 : ! **************************************************************************************************
81 34 : SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
82 : fm_mat_Q_static_bse_gemm, &
83 : Eigenval, Eigenval_scf, &
84 34 : homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
85 34 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
86 :
87 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
88 : fm_mat_S_ab_bse
89 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse_gemm
90 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
91 : INTENT(IN) :: Eigenval, Eigenval_scf
92 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
93 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, bse_lev_virt
94 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
95 : INTEGER, INTENT(IN) :: color_sub
96 : TYPE(mp2_type) :: mp2_env
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
99 : INTEGER, INTENT(IN) :: unit_nr
100 :
101 : CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
102 :
103 : INTEGER :: handle, homo_red, virtual_red
104 : LOGICAL :: my_do_abba, my_do_fulldiag, &
105 : my_do_iterat_diag, my_do_tda
106 : REAL(KIND=dp) :: diag_runtime_est
107 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced
108 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
109 34 : B_bar_ijQ_bse_local, B_iaQ_bse_local
110 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, &
111 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, &
112 : fm_sqrt_A_minus_B
113 : TYPE(cp_logger_type), POINTER :: logger
114 : TYPE(dft_control_type), POINTER :: dft_control
115 : TYPE(mp_para_env_type), POINTER :: para_env
116 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
117 :
118 34 : CALL timeset(routineN, handle)
119 :
120 34 : para_env => fm_mat_S_ia_bse%matrix_struct%para_env
121 :
122 34 : my_do_fulldiag = .FALSE.
123 34 : my_do_iterat_diag = .FALSE.
124 34 : my_do_tda = .FALSE.
125 34 : my_do_abba = .FALSE.
126 : !Method: Iterative or full diagonalization
127 34 : SELECT CASE (mp2_env%bse%bse_diag_method)
128 : CASE (bse_iterdiag)
129 0 : my_do_iterat_diag = .TRUE.
130 : !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
131 0 : CPABORT("Iterative BSE not yet implemented")
132 : CASE (bse_fulldiag)
133 34 : my_do_fulldiag = .TRUE.
134 : END SELECT
135 : !Approximation: TDA and/or full ABBA matrix
136 50 : SELECT CASE (mp2_env%bse%flag_tda)
137 : CASE (bse_tda)
138 16 : my_do_tda = .TRUE.
139 : CASE (bse_abba)
140 18 : my_do_abba = .TRUE.
141 : CASE (bse_both)
142 0 : my_do_tda = .TRUE.
143 34 : my_do_abba = .TRUE.
144 : END SELECT
145 :
146 34 : CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
147 :
148 : ! Link BSE debug flag against debug print level
149 34 : logger => cp_get_default_logger()
150 34 : IF (logger%iter_info%print_level == debug_print_level) THEN
151 0 : mp2_env%bse%bse_debug_print = .TRUE.
152 : END IF
153 :
154 34 : CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
155 : ! We apply the BSE cutoffs using the DFT Eigenenergies
156 : ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
157 : CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
158 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
159 : Eigenval_scf(:, 1, 1), Eigenval(:, 1, 1), Eigenval_reduced, &
160 : homo(1), virtual(1), dimen_RI, unit_nr, &
161 : bse_lev_virt, &
162 : homo_red, virtual_red, &
163 34 : mp2_env)
164 : ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
165 : ! r,s: MO-index, P,R,T: RI-index
166 : ! B: fm_mat_S_..., W: fm_mat_Q_...
167 : CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, &
168 : fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
169 34 : dimen_RI_red, homo_red, virtual_red)
170 :
171 34 : IF (my_do_iterat_diag) THEN
172 : CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, &
173 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
174 : B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
175 : B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, &
176 : gd_array, color_sub, para_env)
177 : END IF
178 :
179 34 : CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env, qs_env)
180 :
181 34 : IF (my_do_fulldiag) THEN
182 : ! Quick estimate of memory consumption and runtime of diagonalizations
183 : CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
184 34 : para_env, diag_runtime_est)
185 : ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
186 : ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
187 : ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
188 : ! α is a spin-dependent factor
189 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
190 : ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
191 :
192 : ! For unscreened W matrix, we need fm_mat_S_ij_trunc
193 34 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
194 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
195 : CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
196 : fm_A_BSE, Eigenval_reduced, unit_nr, &
197 : homo_red, virtual_red, dimen_RI, mp2_env, &
198 4 : para_env, qs_env)
199 : ELSE
200 : CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
201 : fm_A_BSE, Eigenval_reduced, unit_nr, &
202 : homo_red, virtual_red, dimen_RI, mp2_env, &
203 30 : para_env, qs_env)
204 : END IF
205 34 : IF (my_do_abba) THEN
206 : ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
207 : ! B_ia,jb = α * v_ia,jb - W_ib,aj
208 : ! α is a spin-dependent factor
209 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
210 : ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
211 :
212 : ! For unscreened W matrix, we need fm_mat_S_ia_trunc
213 18 : IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
214 : mp2_env%bse%screening_method == bse_screening_alpha) THEN
215 : CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_ia_trunc, fm_B_BSE, &
216 4 : homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
217 : ELSE
218 : CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, &
219 14 : homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
220 : END IF
221 : ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
222 : ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
223 : ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
224 : ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
225 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
226 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
227 18 : homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
228 : END IF
229 34 : CALL cp_fm_release(fm_B_BSE)
230 :
231 34 : NULLIFY (dft_control, tddfpt_control)
232 34 : CALL get_qs_env(qs_env, dft_control=dft_control)
233 34 : tddfpt_control => dft_control%tddfpt2_control
234 34 : IF ((my_do_tda) .AND. (.NOT. tddfpt_control%do_bse)) THEN
235 : ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
236 : CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), &
237 14 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
238 : END IF
239 : ! Release to avoid faulty use of changed A matrix
240 34 : CALL cp_fm_release(fm_A_BSE)
241 34 : IF (my_do_abba) THEN
242 : ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
243 : ! Here, the eigenvectors Z^n relate to X^n via
244 : ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
245 : CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), &
246 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
247 18 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
248 : END IF
249 : ! Release to avoid faulty use of changed C matrix
250 34 : CALL cp_fm_release(fm_C_BSE)
251 : END IF
252 :
253 : CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
254 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
255 34 : fm_mat_Q_static_bse_gemm, mp2_env)
256 34 : DEALLOCATE (Eigenval_reduced)
257 34 : IF (my_do_iterat_diag) THEN
258 : ! Contains untested Block-Davidson algorithm
259 : CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
260 : B_iaQ_bse_local, homo(1), virtual(1), mp2_env%bse%bse_spin_config, unit_nr, &
261 0 : Eigenval(:, 1, 1), para_env, mp2_env)
262 : ! Deallocate local 3c-B-matrices
263 0 : DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
264 : END IF
265 :
266 34 : IF (unit_nr > 0) THEN
267 17 : WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
268 : END IF
269 :
270 34 : CALL timestop(handle)
271 :
272 68 : END SUBROUTINE start_bse_calculation
273 :
274 : END MODULE bse_main
|
