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 Iterative solution of the Bethe-Salpeter equation: the block Davidson solvers on top of
10 : !> the matrix-free application of A and B, and the post-processing of their eigenvectors.
11 : !> The solvers live in bse_davidson, apart from this driver, so that a later Lanczos-Haydock
12 : !> solver for the spectrum can share the matrix-free A and B of bse_matvec and this driver
13 : !> \par History
14 : !> 09.2026 created [Maximilian Graml]
15 : ! **************************************************************************************************
16 : MODULE bse_iterative
17 :
18 : USE bse_davidson, ONLY: abba_indefinite,&
19 : bse_davidson_abba_mk,&
20 : bse_davidson_abba_os,&
21 : bse_davidson_refcheck,&
22 : bse_davidson_tda
23 : USE bse_full_diag, ONLY: create_A_and_B,&
24 : postprocess_bse
25 : USE bse_matvec, ONLY: bse_matvec_create,&
26 : bse_matvec_env_type,&
27 : bse_matvec_release,&
28 : bse_matvec_selfcheck
29 : USE bse_print, ONLY: print_output_header
30 : USE bse_types, ONLY: bse_env_type
31 : USE cp_blacs_env, ONLY: cp_blacs_env_release,&
32 : cp_blacs_env_type
33 : USE cp_control_types, ONLY: dft_control_type
34 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 : cp_fm_struct_release,&
36 : cp_fm_struct_type
37 : USE cp_fm_types, ONLY: cp_fm_create,&
38 : cp_fm_get_info,&
39 : cp_fm_release,&
40 : cp_fm_to_fm_submat_general,&
41 : cp_fm_type
42 : USE input_constants, ONLY: bse_abba_mk_davidson,&
43 : bse_abba_olsen_stratmann,&
44 : bse_screening_alpha,&
45 : bse_screening_rpa,&
46 : bse_screening_tdhf,&
47 : bse_singlet,&
48 : bse_triplet
49 : USE kinds, ONLY: dp
50 : USE physcon, ONLY: evolt
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_iterative'
58 :
59 : PUBLIC :: solve_bse_iteratively
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief Lowest excitations of a closed-shell reference from the block Davidson solvers on top of
65 : !> the matrix-free application of A and B, post-processed as in the full diagonalization.
66 : !> TDA (bse_davidson_tda): sum_jb A_ia,jb X_jb^n = Ω^n X_ia^n. Full problem, by ABBA_SOLVER:
67 : !> MK_DAVIDSON, sum_jb [(A+B)(A-B)]_ia,jb x_jb^n = (Ω^n)^2 x_ia^n with x^n = X^n - Y^n;
68 : !> OLSEN_STRATMANN, sum_jb (A+B)_ia,jb y_jb^n = Ω^n x_ia^n and sum_jb (A-B)_ia,jb x_jb^n = Ω^n y_ia^n
69 : !> with y^n = X^n + Y^n. Neither A nor B is formed: bse_matvec_apply contracts the RI slabs,
70 : !> (A Z)_ia = (ε_a-ε_i) Z_ia + α sum_P B^P_ia sum_jb B^P_jb Z_jb - w sum_Pjb \bar{B}^P_ij B^P_ab Z_jb
71 : !> (B Z)_ia = α sum_P B^P_ia sum_jb B^P_jb Z_jb - w sum_Pjb \bar{B}^P_ib B^P_ja Z_jb
72 : !> with α, w and \bar{B} set by create_matvec_env from SPIN_CONFIG and SCREENING_IN_W.
73 : !> The five slabs are released once the solver holds its RI-sliced copies
74 : !> \param Eigenval_reduced energies of the homo_red + virt_red active levels in Hartree, QP or DFT per USE_KS_ENERGIES
75 : !> \param homo_red_arr occupied levels in the active window, per spin; entry 1 is used
76 : !> \param virt_red_arr virtual levels in the active window, per spin; entry 1 is used
77 : !> \param homo occupied levels of the full reference, per spin; postprocess_bse indexes the MOs with it
78 : !> \param dimen_RI number of RI functions; only the debug construction of the explicit A and B uses it
79 : !> \param do_tda solve the TDA problem
80 : !> \param do_abba solve the full problem; after the TDA in a TDA+ABBA run, seeded with X^TDA
81 : !> \param bse_env the BSE environment: the BSE_ITERAT settings, BSE_DEBUG_PRINT and the TDDFPT control flags;
82 : !> the five slabs of spin 1, released here: fm_S_ia = B^P_ia (N_RI x homo_red*virt_red), fm_S_ij = B^P_ij,
83 : !> fm_S_ab = B^P_ab, fm_S_bar_ia = \bar{B}^P_ia = sum_Q [1+Q(0)]^-1_PQ B^Q_ia (the slab of W_ib,aj) and
84 : !> fm_S_bar_ij = \bar{B}^P_ij (the slab of W_ij,ab)
85 : !> \param mo_coeff MO coefficients of the reference, handed to postprocess_bse
86 : !> \param unit_nr output unit, positive on the writing rank only
87 : ! **************************************************************************************************
88 24 : SUBROUTINE solve_bse_iteratively(Eigenval_reduced, &
89 8 : homo_red_arr, virt_red_arr, homo, dimen_RI, &
90 8 : do_tda, do_abba, bse_env, mo_coeff, unit_nr)
91 :
92 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval_reduced
93 : INTEGER, DIMENSION(:), INTENT(IN) :: homo_red_arr, virt_red_arr, homo
94 : INTEGER, INTENT(IN) :: dimen_RI
95 : LOGICAL, INTENT(IN) :: do_tda, do_abba
96 : TYPE(bse_env_type), INTENT(INOUT) :: bse_env
97 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
98 : INTEGER, INTENT(IN) :: unit_nr
99 :
100 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_bse_iteratively'
101 : INTEGER, PARAMETER :: n_ov_max_debug = 20000
102 :
103 : CHARACTER(LEN=16) :: margin_str, n_ov_str, solver_str
104 : INTEGER :: abba_status, handle, n_ov
105 : LOGICAL :: do_debug
106 : REAL(KIND=dp) :: ab_margin
107 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens, exc_ens_abba
108 8 : TYPE(bse_matvec_env_type) :: mv_env
109 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_ia
110 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_X_abba, fm_X_ia, &
111 : fm_Y_abba, fm_Y_ia
112 : TYPE(cp_fm_type), POINTER :: fm_tda_guess
113 : TYPE(cp_fm_type), TARGET :: fm_X
114 : TYPE(dft_control_type), POINTER :: dft_control
115 :
116 8 : CALL timeset(routineN, handle)
117 :
118 8 : dft_control => bse_env%dft_control
119 8 : IF (dft_control%tddfpt2_control%do_bse .OR. dft_control%tddfpt2_control%do_bse_w_only .OR. &
120 : dft_control%tddfpt2_control%do_bse_gw_only) THEN
121 : ! the TDDFPT path reads the explicit A and W matrices that create_A_and_B (FULLDIAG) stashes in the
122 : ! excited-state environment; ITERDIAG never forms them
123 : CALL cp_abort(__LOCATION__, &
124 : "BSE_DIAG_METHOD ITERDIAG cannot provide the BSE kernel of TDDFPT "// &
125 0 : "(DO_BSE, DO_BSE_W_ONLY, DO_BSE_GW_ONLY). Use BSE_DIAG_METHOD FULLDIAG.")
126 : END IF
127 :
128 8 : n_ov = homo_red_arr(1)*virt_red_arr(1)
129 : ! the debug checks hold the explicit N_ov x N_ov matrices
130 8 : do_debug = bse_env%bse_debug_print .AND. n_ov <= n_ov_max_debug
131 8 : IF (bse_env%bse_debug_print .AND. .NOT. do_debug .AND. unit_nr > 0) THEN
132 0 : WRITE (n_ov_str, '(I0)') n_ov_max_debug
133 : CALL cp_warn(__LOCATION__, &
134 : "BSE_DEBUG_PRINT: the checks of the iterative solver against the explicit A "// &
135 0 : "and B are skipped above N_ov = "//TRIM(n_ov_str)//".")
136 : END IF
137 :
138 : CALL create_matvec_env(mv_env, bse_env%fm_S_ia(1), bse_env%fm_S_ij(1), &
139 : bse_env%fm_S_ab(1), bse_env%fm_S_bar_ia(1), &
140 : bse_env%fm_S_bar_ij(1), Eigenval_reduced, &
141 8 : homo_red_arr(1), virt_red_arr(1), do_abba, bse_env, unit_nr)
142 :
143 8 : IF (do_debug) THEN
144 : CALL create_A_and_B(bse_env%fm_S_ia, bse_env%fm_S_ij, bse_env%fm_S_ab, &
145 : bse_env%fm_S_bar_ia, bse_env%fm_S_bar_ij, fm_A_BSE, fm_B_BSE, &
146 : do_abba, Eigenval_reduced, unit_nr, homo_red_arr, virt_red_arr, dimen_RI, &
147 0 : bse_env)
148 0 : IF (do_abba) THEN
149 0 : CALL bse_matvec_selfcheck(mv_env, fm_A_BSE, unit_nr, fm_B_BSE)
150 : ELSE
151 0 : CALL bse_matvec_selfcheck(mv_env, fm_A_BSE, unit_nr)
152 : END IF
153 : END IF
154 :
155 : ! the solver works on its RI-sliced copies; only the grid of the slabs is kept, for the eigenvectors
156 8 : blacs_env_ia => bse_env%fm_S_ia(1)%matrix_struct%context
157 8 : CALL blacs_env_ia%retain()
158 8 : CALL cp_fm_release(bse_env%fm_S_ia(1))
159 8 : CALL cp_fm_release(bse_env%fm_S_ij(1))
160 8 : CALL cp_fm_release(bse_env%fm_S_ab(1))
161 8 : CALL cp_fm_release(bse_env%fm_S_bar_ia(1))
162 8 : CALL cp_fm_release(bse_env%fm_S_bar_ij(1))
163 :
164 8 : IF (do_tda) THEN
165 4 : CALL print_output_header(homo_red_arr(1), virt_red_arr(1), homo(1), .TRUE., bse_env, unit_nr)
166 4 : CALL bse_davidson_tda(mv_env, bse_env, unit_nr, exc_ens, fm_X)
167 4 : IF (do_debug) CALL bse_davidson_refcheck(fm_A_BSE, exc_ens, fm_X, bse_env, unit_nr)
168 :
169 4 : CALL copy_to_ia_context(fm_X, blacs_env_ia, fm_X_ia)
170 : CALL postprocess_bse(exc_ens, fm_X_ia, bse_env, mo_coeff, &
171 4 : homo_red_arr(1), virt_red_arr(1), homo(1), unit_nr, .TRUE.)
172 4 : CALL cp_fm_release(fm_X_ia)
173 4 : DEALLOCATE (exc_ens)
174 : END IF
175 :
176 8 : IF (do_abba) THEN
177 : ! the TDA vectors are the initial guess of a combined run; a null pointer passed to the
178 : ! optional argument counts as absent (Fortran 2008)
179 4 : NULLIFY (fm_tda_guess)
180 4 : IF (do_tda) fm_tda_guess => fm_X
181 4 : CALL print_output_header(homo_red_arr(1), virt_red_arr(1), homo(1), .FALSE., bse_env, unit_nr)
182 6 : SELECT CASE (bse_env%abba_solver)
183 : CASE (bse_abba_mk_davidson)
184 2 : solver_str = 'MK_DAVIDSON'
185 : CALL bse_davidson_abba_mk(mv_env, bse_env, unit_nr, exc_ens_abba, fm_X_abba, fm_Y_abba, &
186 2 : abba_status, ab_margin, fm_tda_guess)
187 : CASE (bse_abba_olsen_stratmann)
188 2 : solver_str = 'OLSEN_STRATMANN'
189 : CALL bse_davidson_abba_os(mv_env, bse_env, unit_nr, exc_ens_abba, fm_X_abba, fm_Y_abba, &
190 2 : abba_status, ab_margin, fm_tda_guess)
191 : CASE DEFAULT
192 4 : CPABORT("Unknown BSE_ITERAT%ABBA_SOLVER")
193 : END SELECT
194 4 : IF (abba_status == abba_indefinite) THEN
195 0 : WRITE (margin_str, '(F16.6)') ab_margin*evolt
196 0 : IF (.NOT. do_tda) THEN
197 : CALL cp_abort(__LOCATION__, &
198 : "Matrix (A-B) is not positive definite on the trial space of the "// &
199 : TRIM(solver_str)//" solver (eigenvalue "//TRIM(ADJUSTL(margin_str))//" eV). "// &
200 0 : "The full ABBA problem is ill-defined for this reference.")
201 : END IF
202 : ! the TDA results of a combined run stand on their own
203 0 : IF (unit_nr > 0) THEN
204 : CALL cp_warn(__LOCATION__, &
205 : "BSE: A-B is not positive definite on the trial space (margin "// &
206 : TRIM(ADJUSTL(margin_str))//" eV); the ABBA calculation is skipped, "// &
207 0 : "TDA results above are complete.")
208 0 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
209 0 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.6)') 'BSE|', &
210 0 : 'A-B is not positive definite on the trial space (eV)', ab_margin*evolt
211 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
212 0 : 'The ABBA calculation is skipped, the TDA results above are complete'
213 : END IF
214 : ELSE
215 4 : IF (do_debug) THEN
216 : CALL bse_davidson_refcheck(fm_A_BSE, exc_ens_abba, fm_X_abba, bse_env, unit_nr, &
217 0 : fm_B_BSE, fm_Y_abba)
218 : END IF
219 :
220 4 : CALL copy_to_ia_context(fm_X_abba, blacs_env_ia, fm_X_ia)
221 4 : CALL copy_to_ia_context(fm_Y_abba, blacs_env_ia, fm_Y_ia)
222 4 : CALL cp_fm_release(fm_X_abba)
223 4 : CALL cp_fm_release(fm_Y_abba)
224 : CALL postprocess_bse(exc_ens_abba, fm_X_ia, bse_env, mo_coeff, &
225 4 : homo_red_arr(1), virt_red_arr(1), homo(1), unit_nr, .FALSE., fm_Y_ia)
226 4 : CALL cp_fm_release(fm_X_ia)
227 4 : CALL cp_fm_release(fm_Y_ia)
228 4 : DEALLOCATE (exc_ens_abba)
229 : END IF
230 : END IF
231 :
232 8 : IF (do_tda) CALL cp_fm_release(fm_X)
233 8 : IF (do_debug) THEN
234 0 : CALL cp_fm_release(fm_A_BSE)
235 0 : IF (do_abba) CALL cp_fm_release(fm_B_BSE)
236 : END IF
237 8 : CALL cp_blacs_env_release(blacs_env_ia)
238 8 : CALL bse_matvec_release(mv_env)
239 :
240 8 : CALL timestop(handle)
241 :
242 16 : END SUBROUTINE solve_bse_iteratively
243 :
244 : ! **************************************************************************************************
245 : !> \brief Copies eigenvectors from the process grid of the solver to the one of the RI slabs, from
246 : !> which postprocess_bse takes its grid
247 : !> \param fm_vec eigenvectors on the npe x 1 process grid of the solver
248 : !> \param blacs_env_ia process grid of the RI slabs
249 : !> \param fm_vec_ia the copy on blacs_env_ia, created here
250 : ! **************************************************************************************************
251 24 : SUBROUTINE copy_to_ia_context(fm_vec, blacs_env_ia, fm_vec_ia)
252 :
253 : TYPE(cp_fm_type), INTENT(IN) :: fm_vec
254 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_ia
255 : TYPE(cp_fm_type), INTENT(OUT) :: fm_vec_ia
256 :
257 : INTEGER :: n_exc, n_ov
258 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
259 :
260 12 : CALL cp_fm_get_info(fm_vec, nrow_global=n_ov, ncol_global=n_exc)
261 12 : NULLIFY (fm_struct)
262 : CALL cp_fm_struct_create(fm_struct, para_env=blacs_env_ia%para_env, context=blacs_env_ia, &
263 12 : nrow_global=n_ov, ncol_global=n_exc)
264 12 : CALL cp_fm_create(fm_vec_ia, fm_struct, name="fm_davidson_ia")
265 12 : CALL cp_fm_struct_release(fm_struct)
266 12 : CALL cp_fm_to_fm_submat_general(fm_vec, fm_vec_ia, n_ov, n_exc, 1, 1, 1, 1, blacs_env_ia)
267 :
268 12 : END SUBROUTINE copy_to_ia_context
269 :
270 : ! **************************************************************************************************
271 : !> \brief Sets up the matrix-free application of A and B with the prefactors and the screened
272 : !> slabs that create_A and create_B use for the chosen SPIN_CONFIG and SCREENING_IN_W:
273 : !> α = 2 (SINGLET) or 0 (TRIPLET); w = 1 with the W-screened \bar{B} (default screening),
274 : !> w = SCREENING_FACTOR with \bar{B} = B (ALPHA), w = 1 with \bar{B} = B (TDHF), w = 0 (RPA)
275 : !> \param mv_env the environment created here
276 : !> \param fm_S_ia B^P_ia, N_RI x homo_red*virt_red
277 : !> \param fm_S_ij B^P_ij, N_RI x homo_red^2
278 : !> \param fm_S_ab B^P_ab, N_RI x virt_red^2
279 : !> \param fm_S_bar_ia W-screened \bar{B}^P_ia
280 : !> \param fm_S_bar_ij W-screened \bar{B}^P_ij
281 : !> \param Eigenval_reduced energies of the homo_red + virt_red active levels in Hartree, QP or DFT per USE_KS_ENERGIES
282 : !> \param homo_red occupied levels in the active window
283 : !> \param virt_red virtual levels in the active window
284 : !> \param do_abba B is applied too, so the \bar{B}^P_ia slab is sliced as well
285 : !> \param bse_env SPIN_CONFIG, SCREENING_IN_W, SCREENING_FACTOR and MATVEC_BLOCK_COLUMNS
286 : !> \param unit_nr output unit, positive on the writing rank only
287 : ! **************************************************************************************************
288 8 : SUBROUTINE create_matvec_env(mv_env, fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, fm_S_bar_ij, &
289 8 : Eigenval_reduced, homo_red, virt_red, do_abba, bse_env, unit_nr)
290 :
291 : TYPE(bse_matvec_env_type), INTENT(OUT) :: mv_env
292 : TYPE(cp_fm_type), INTENT(IN), TARGET :: fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, &
293 : fm_S_bar_ij
294 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval_reduced
295 : INTEGER, INTENT(IN) :: homo_red, virt_red
296 : LOGICAL, INTENT(IN) :: do_abba
297 : TYPE(bse_env_type), INTENT(IN) :: bse_env
298 : INTEGER, INTENT(IN) :: unit_nr
299 :
300 : REAL(KIND=dp) :: alpha, w_fac
301 : TYPE(cp_fm_type), POINTER :: fm_W_ia, fm_W_ij
302 :
303 8 : alpha = 0.0_dp
304 16 : SELECT CASE (bse_env%bse_spin_config)
305 : CASE (bse_singlet)
306 8 : alpha = 2.0_dp
307 : CASE (bse_triplet)
308 8 : alpha = 0.0_dp
309 : END SELECT
310 :
311 8 : SELECT CASE (bse_env%screening_method)
312 : CASE (bse_screening_alpha)
313 0 : w_fac = bse_env%screening_factor
314 : CASE (bse_screening_rpa)
315 0 : w_fac = 0.0_dp
316 : CASE DEFAULT
317 8 : w_fac = 1.0_dp
318 : END SELECT
319 :
320 : ! TDHF and ALPHA screening contract the bare slabs, all others the W-multiplied ones
321 8 : IF (bse_env%screening_method == bse_screening_tdhf .OR. &
322 : bse_env%screening_method == bse_screening_alpha) THEN
323 : fm_W_ij => fm_S_ij
324 : fm_W_ia => fm_S_ia
325 : ELSE
326 8 : fm_W_ij => fm_S_bar_ij
327 8 : fm_W_ia => fm_S_bar_ia
328 : END IF
329 : CALL bse_matvec_create(mv_env, fm_S_ia, fm_W_ij, fm_S_ab, fm_W_ia, Eigenval_reduced, &
330 : homo_red, virt_red, alpha, w_fac, do_abba, unit_nr, &
331 8 : block_cols=bse_env%matvec_block_columns, memory_per_proc_gb=bse_env%memory_per_proc_gb)
332 :
333 8 : END SUBROUTINE create_matvec_env
334 :
335 : END MODULE bse_iterative
|
