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 Block Davidson solvers for the lowest excitations of the Bethe-Salpeter equation on top
10 : !> of the matrix-free application of A and B: the TDA problem sum_jb A_ia,jb X_jb = Ω X_ia
11 : !> (bse_davidson_tda) and the full problem, either as sum_jb [(A+B)(A-B)]_ia,jb x_jb = Ω^2 x_ia
12 : !> (bse_davidson_abba_mk) or as the pair sum_jb (A+B)_ia,jb y_jb = Ω x_ia, sum_jb (A-B)_ia,jb x_jb = Ω y_ia
13 : !> (bse_davidson_abba_os). Every solver projects its problem onto an orthonormal basis, takes the
14 : !> lowest Ritz pairs, extends the basis by the preconditioned residuals of the roots still in need
15 : !> and restarts thick at the subspace ceiling; the helpers below are shared by the three drivers
16 : !> \par History
17 : !> 09.2026 created [Maximilian Graml]
18 : ! **************************************************************************************************
19 : MODULE bse_davidson
20 :
21 : USE bibliography, ONLY: Bai2012,&
22 : Davidson1975,&
23 : Fukaya2014,&
24 : Olsen1988,&
25 : Stathopoulos1998,&
26 : Stratmann1998,&
27 : Vecharynski2017,&
28 : cite_reference
29 : USE bse_full_diag, ONLY: create_hermitian_form_of_ABBA
30 : USE bse_matvec, ONLY: bse_matvec_apply,&
31 : bse_matvec_diagonal,&
32 : bse_matvec_env_type,&
33 : bse_matvec_subblock,&
34 : bse_matvec_vector_struct,&
35 : mem_fraction
36 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
37 : cp_blacs_env_release,&
38 : cp_blacs_env_type
39 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
40 : cp_fm_scale_and_add,&
41 : cp_fm_transpose
42 : USE cp_fm_diag, ONLY: choose_eigv_solver
43 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
44 : cp_fm_struct_release,&
45 : cp_fm_struct_type
46 : USE cp_fm_types, ONLY: &
47 : cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
48 : cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type
49 : USE input_constants, ONLY: bse_iter_and_cond,&
50 : bse_iter_en_cond,&
51 : bse_iter_or_cond,&
52 : bse_iter_res_cond,&
53 : bse_memcheck_abort,&
54 : bse_memcheck_clamp,&
55 : bse_memcheck_off,&
56 : bse_memcheck_warn
57 : USE kinds, ONLY: dp
58 : USE machine, ONLY: m_walltime
59 : USE mathlib, ONLY: diamat_all
60 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
61 : mp_para_env_type
62 : USE mp2_types, ONLY: mp2_type
63 : USE parallel_gemm_api, ONLY: parallel_gemm
64 : USE physcon, ONLY: evolt
65 : USE util, ONLY: sort
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_davidson'
73 :
74 : ! energies closer than this belong to one multiplet (Hartree)
75 : REAL(KIND=dp), PARAMETER, PRIVATE :: deg_thresh = 1.0E-6_dp
76 : ! a basis that is not orthonormal to this accuracy makes the residuals meaningless
77 : REAL(KIND=dp), PARAMETER, PRIVATE :: ortho_tol = 1.0E-8_dp
78 :
79 : ! outcome of an ABBA solver: solved, or A-B not positive definite on the trial space
80 : INTEGER, PARAMETER, PUBLIC :: abba_ok = 0, abba_indefinite = 1
81 :
82 : ! the three drivers, for the memory estimate
83 : INTEGER, PARAMETER, PRIVATE :: driver_tda = 1, driver_mk = 2, driver_os = 3
84 :
85 : PUBLIC :: bse_davidson_tda, bse_davidson_abba_mk, bse_davidson_abba_os, &
86 : bse_davidson_refcheck
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of
92 : !> sum_jb A_ia,jb X_jb^n = Ω^n X_ia^n. With the orthonormal basis Z_ia,m the reduced matrix is
93 : !> Ar_mn = sum_ia Z_ia,m (A Z)_ia,n with eigenpairs (θ_k, c_mk), X_ia^k = sum_m Z_ia,m c_mk,
94 : !> and the basis grows by the corrections t_ia,k = r_ia,k/(d_ia - θ_k) of the residuals
95 : !> r_ia,k = sum_m (A Z)_ia,m c_mk - θ_k X_ia^k, d_ia being the diagonal chosen by PRECONDITIONER.
96 : !> Arrays: Ar_mn (the reduced A) in fm_red_A on the process grid, c_mk in coef_ritz, θ_k c_mk in
97 : !> coef_ritz_theta
98 : !> \param mv_env slabs, prefactors and transition energies of the matrix-free A
99 : !> \param mp2_env mp2_env%bse carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
100 : !> \param unit_nr output unit, positive on the writing rank only
101 : !> \param exc_ens Ω^n in Hartree
102 : !> \param fm_X X_ia^n on the process grid of mv_env, column n
103 : ! **************************************************************************************************
104 2 : SUBROUTINE bse_davidson_tda(mv_env, mp2_env, unit_nr, exc_ens, fm_X)
105 :
106 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
107 : TYPE(mp2_type), INTENT(IN) :: mp2_env
108 : INTEGER, INTENT(IN) :: unit_nr
109 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
110 : INTENT(OUT) :: exc_ens
111 : TYPE(cp_fm_type), INTENT(OUT) :: fm_X
112 :
113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_davidson_tda'
114 :
115 : INTEGER :: block_size, handle, iter, k, m, m_max, &
116 : n_act, n_cand, n_dependent, n_kernel, &
117 : n_ov, n_req, n_restart, n_want, nt
118 : INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
119 : LOGICAL :: all_conv, has_prev, stalled
120 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
121 : REAL(KIND=dp) :: dev, t_iter
122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, diag, res, theta, theta_prev, &
123 2 : theta_sub
124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_restart, coef_restart_cand, &
125 2 : coef_ritz, coef_ritz_theta, red_block
126 : TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
127 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
128 : TYPE(cp_fm_type) :: fm_AZ, fm_red_A, fm_work, fm_Z
129 : TYPE(mp_para_env_type), POINTER :: para_env
130 :
131 2 : CALL timeset(routineN, handle)
132 2 : CALL cite_reference(Davidson1975)
133 2 : CALL cite_reference(Stathopoulos1998)
134 2 : CALL cite_reference(Fukaya2014)
135 :
136 2 : para_env => mv_env%para_env
137 : ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
138 2 : CPASSERT(mv_env%blacs_env%num_pe(2) == 1)
139 : ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
140 2 : NULLIFY (diag_blacs_env)
141 2 : CALL cp_blacs_env_create(diag_blacs_env, para_env)
142 2 : n_ov = mv_env%n_ov
143 2 : CALL davidson_sizes(mp2_env, n_ov, unit_nr, n_want, n_act, block_size)
144 :
145 6 : ALLOCATE (diag(n_ov))
146 2 : CALL bse_matvec_diagonal(mv_env, mp2_env%bse%preconditioner, diag)
147 :
148 2 : CALL subspace_ceiling(mp2_env, mv_env, driver_tda, n_act, block_size, unit_nr, m_max)
149 :
150 2 : CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
151 2 : CALL cp_fm_create(fm_Z, fm_struct, name="fm_Z_davidson", set_zero=.TRUE.)
152 2 : CALL cp_fm_create(fm_AZ, fm_struct, name="fm_AZ_davidson", set_zero=.TRUE.)
153 2 : CALL cp_fm_struct_release(fm_struct)
154 2 : CALL bse_matvec_vector_struct(mv_env, MIN(2*n_act, m_max), fm_struct)
155 2 : CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.TRUE.)
156 2 : CALL cp_fm_struct_release(fm_struct)
157 2 : NULLIFY (fm_struct)
158 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
159 2 : nrow_global=m_max, ncol_global=m_max)
160 2 : CALL cp_fm_create(fm_red_A, fm_struct, name="fm_red_A_davidson", set_zero=.TRUE.)
161 2 : CALL cp_fm_struct_release(fm_struct)
162 :
163 : ! c_mk and θ_k c_mk, the restart candidates (current and previous Ritz vectors) and the restart basis Q_pk
164 20 : ALLOCATE (coef_ritz(m_max, n_act), coef_restart_cand(m_max, 2*n_act), coef_ritz_theta(m_max, n_act))
165 6 : ALLOCATE (coef_restart(m_max, 2*n_act))
166 : ALLOCATE (theta(m_max), theta_prev(n_act), res(n_act), dE(n_act), conv(n_act), need(n_act), &
167 22 : selected(n_act))
168 2 : coef_restart_cand(:, :) = 0.0_dp
169 2 : theta_prev(:) = 0.0_dp
170 2 : has_prev = .FALSE.
171 2 : stalled = .FALSE.
172 2 : n_restart = 0
173 2 : n_dependent = 0
174 :
175 2 : IF (mp2_env%bse%num_guess_transitions > 0) THEN
176 : CALL initial_guess_subblock(mv_env, mp2_env, diag, n_act, m_max, deg_thresh, fm_Z, m, theta_sub, &
177 0 : unit_nr)
178 : ELSE
179 2 : CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_Z, m)
180 : END IF
181 : ! both guess routines return orthonormal columns by construction; one Gram matrix checks it
182 2 : IF (extension_deviation(fm_Z, fm_Z, 0, m, para_env) > ortho_tol) THEN
183 0 : CPABORT("BSE Davidson: the initial guess is not orthonormal")
184 : END IF
185 : ! (A Z)_ia,m of the guess and Ar_mn = sum_ia Z_ia,m (A Z)_ia,n
186 2 : CALL bse_matvec_apply(mv_env, fm_Z, 1, m, fm_AZ)
187 2 : n_kernel = m
188 2 : CALL reduced_gram_blocks(fm_Z, m, fm_AZ, 1, m, block_size, para_env, fm_red_A)
189 :
190 2 : CALL print_iteration_header('Block Davidson iterations within the TDA:', unit_nr)
191 :
192 2 : all_conv = .FALSE.
193 16 : DO iter = 1, mp2_env%bse%max_iter
194 16 : t_iter = m_walltime()
195 :
196 : ! sum_n Ar_mn c_nk = θ_k c_mk, θ ascending
197 16 : CALL solve_reduced(fm_red_A, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
198 :
199 : ! thick restart onto the current and the previous Ritz vectors, no kernel application; with the
200 : ! restart basis Q_pk in coef_restart, Z_ia,k <- sum_p Z_ia,p Q_pk, A Z alike, Ar_kl <- sum_pq Q_pk Ar_pq Q_ql
201 16 : IF (m + block_size > m_max .AND. m > n_act .AND. m < n_ov) THEN
202 706 : coef_restart_cand(1:m, 1:n_act) = coef_ritz(1:m, 1:n_act)
203 6 : n_cand = n_act
204 6 : IF (has_prev) n_cand = 2*n_act
205 6 : CALL restart_basis(coef_restart_cand, m, n_cand, 2*n_act, coef_restart, k)
206 6 : CALL rotate_in_place(fm_Z, m, coef_restart, k, fm_work, fm_AZ)
207 6 : CALL reduced_rotate(fm_red_A, m, coef_restart, k, para_env, diag_blacs_env)
208 6 : m = k
209 6 : n_restart = n_restart + 1
210 6 : CALL solve_reduced(fm_red_A, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
211 : END IF
212 :
213 : ! r_ia,k = sum_m (A Z)_ia,m c_mk - θ_k sum_m Z_ia,m c_mk for all tracked roots k
214 128 : DO k = 1, n_act
215 1444 : coef_ritz_theta(1:m, k) = coef_ritz(1:m, k)*theta(k)
216 : END DO
217 1444 : CALL subspace_rotate(fm_AZ, m, coef_ritz(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
218 1444 : CALL subspace_rotate(fm_Z, m, coef_ritz_theta(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
219 16 : CALL column_norms(fm_work, 1, n_act, para_env, res)
220 :
221 : CALL convergence_step(mp2_env, theta, theta_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
222 16 : dE, conv, n_req, need)
223 42 : all_conv = .NOT. ANY(need)
224 16 : IF (all_conv) EXIT
225 :
226 14 : IF (iter == mp2_env%bse%max_iter) EXIT
227 :
228 : ! corrections for the lowest unconverged roots, t_ia,k = r_ia,k/(d_ia - θ_k)
229 14 : CALL select_roots(mp2_env, need, res, MIN(block_size, m_max - m), selected, nt)
230 14 : CALL davidson_corrections(fm_work, selected(1:nt), theta, diag, fm_Z, m + 1)
231 14 : CALL extend_orthonormal(fm_Z, m, nt, para_env, n_dependent)
232 : ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
233 : ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
234 14 : IF (nt == 0 .AND. stalled) THEN
235 0 : CPABORT("BSE Davidson: correction vectors vanished before convergence")
236 : END IF
237 14 : stalled = nt == 0
238 : ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
239 14 : IF (extension_deviation(fm_Z, fm_Z, m, nt, para_env) > ortho_tol) THEN
240 0 : CPABORT("BSE Davidson: the basis lost orthonormality")
241 : END IF
242 :
243 14 : IF (nt > 0) THEN
244 14 : CALL bse_matvec_apply(mv_env, fm_Z, m + 1, nt, fm_AZ)
245 14 : n_kernel = n_kernel + nt
246 :
247 : ! new columns Ar_mn = sum_ia Z_ia,m (A Z)_ia,n of the symmetric reduced matrix
248 56 : ALLOCATE (red_block(m + nt, nt))
249 14 : CALL subspace_gram(fm_Z, 1, m + nt, fm_AZ, m + 1, nt, para_env, red_block)
250 14 : CALL reduced_extend(fm_red_A, red_block, m, nt)
251 14 : DEALLOCATE (red_block)
252 : END IF
253 :
254 : ! the Ritz vectors of this iteration are the second block of restart candidates
255 2170 : coef_restart_cand(:, n_act + 1:2*n_act) = 0.0_dp
256 1260 : coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_ritz(1:m, 1:n_act)
257 112 : has_prev = .TRUE.
258 112 : theta_prev(:) = theta(1:n_act)
259 32 : m = m + nt
260 : END DO
261 :
262 2 : IF (.NOT. all_conv) CALL abort_unconverged(theta, res, conv, need, n_req, unit_nr)
263 2 : IF (mp2_env%bse%bse_debug_print) CALL print_tracked_roots(theta, res, n_want, n_act, unit_nr)
264 :
265 2 : IF (mp2_env%bse%bse_debug_print) THEN
266 0 : dev = orthonormality_deviation(fm_Z, fm_Z, m, block_size, para_env)
267 0 : IF (unit_nr > 0) THEN
268 0 : WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
269 0 : 'Max deviation of the basis from orthonormality', dev
270 : END IF
271 : END IF
272 :
273 2 : IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, theta, n_req, unit_nr)
274 2 : CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, n_dependent=n_dependent)
275 :
276 6 : ALLOCATE (exc_ens(n_want))
277 8 : exc_ens(:) = theta(1:n_want)
278 2 : CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
279 2 : CALL cp_fm_create(fm_X, fm_struct, name="fm_X_davidson", set_zero=.TRUE.)
280 2 : CALL cp_fm_struct_release(fm_struct)
281 : ! X_ia^n = sum_m Z_ia,m c_mn
282 80 : CALL subspace_rotate(fm_Z, m, coef_ritz(1:m, 1:n_want), n_want, fm_X, 1.0_dp, 0.0_dp)
283 :
284 2 : CALL cp_fm_release(fm_Z)
285 2 : CALL cp_fm_release(fm_AZ)
286 2 : CALL cp_fm_release(fm_work)
287 2 : CALL cp_fm_release(fm_red_A)
288 2 : CALL cp_blacs_env_release(diag_blacs_env)
289 0 : DEALLOCATE (coef_ritz, coef_restart_cand, coef_ritz_theta, coef_restart, theta, theta_prev, res, dE, conv, &
290 2 : need, selected, diag)
291 2 : IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
292 :
293 2 : CALL timestop(handle)
294 :
295 6 : END SUBROUTINE bse_davidson_tda
296 :
297 : ! **************************************************************************************************
298 : !> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem,
299 : !> written as sum_jb [(A+B)(A-B)]_ia,jb x_jb^n = (Ω^n)^2 x_ia^n with x^n = X^n - Y^n, which is
300 : !> symmetric in the inner product <u,v> = sum_ia,jb u_ia (A-B)_ia,jb v_jb. The basis V_ia,m is
301 : !> orthonormal in that inner product, W_ia,m = sum_jb (A-B)_ia,jb V_jb,m, and the reduced matrix
302 : !> Hr_mn = sum_ia,jb W_ia,m (A+B)_ia,jb W_jb,n has the eigenpairs ((Ω^k)^2, c_mk). Then
303 : !> x_ia^k = (Ω^k)^1/2 sum_m V_ia,m c_mk and y_ia^k = (X^k+Y^k)_ia = (Ω^k)^-1/2 sum_m W_ia,m c_mk
304 : !> obey sum_ia x_ia^k y_ia^k = 1, the residual is r_ia,k = sum_jb (A+B)_ia,jb y_jb^k - Ω^k x_ia^k
305 : !> and the correction t_ia,k = r_ia,k/(d_ia^2 - (Ω^k)^2), d_ia chosen by PRECONDITIONER.
306 : !> Arrays: Hr_mn in fm_red_M on the process grid, c_mk in coef_ritz, (Ω^k)^1/2 c_mk in coef_x,
307 : !> (Ω^k)^-1/2 c_mk in coef_y, Ω^k (Ω^k)^1/2 c_mk in coef_x_omega
308 : !> \param mv_env slabs, prefactors and transition energies of the matrix-free A and B
309 : !> \param mp2_env mp2_env%bse carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
310 : !> \param unit_nr output unit, positive on the writing rank only
311 : !> \param exc_ens Ω^n in Hartree
312 : !> \param fm_X X_ia^n on the process grid of mv_env, column n; created only for abba_ok
313 : !> \param fm_Y Y_ia^n, as fm_X
314 : !> \param abba_status abba_indefinite if A-B is not positive definite on the trial space
315 : !> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
316 : !> \param fm_X_tda converged TDA vectors as initial guess
317 : ! **************************************************************************************************
318 2 : SUBROUTINE bse_davidson_abba_mk(mv_env, mp2_env, unit_nr, exc_ens, fm_X, fm_Y, abba_status, &
319 : ab_margin, fm_X_tda)
320 :
321 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
322 : TYPE(mp2_type), INTENT(IN) :: mp2_env
323 : INTEGER, INTENT(IN) :: unit_nr
324 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
325 : INTENT(OUT) :: exc_ens
326 : TYPE(cp_fm_type), INTENT(OUT) :: fm_X, fm_Y
327 : INTEGER, INTENT(OUT) :: abba_status
328 : REAL(KIND=dp), INTENT(OUT) :: ab_margin
329 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_X_tda
330 :
331 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_davidson_abba_mk'
332 :
333 : INTEGER :: block_size, guess_col, handle, iter, k, m, m_max, n_act, n_cand, n_dependent, &
334 : n_guess_cols, n_kernel, n_ov, n_req, n_restart, n_seed, n_want, nrow_local, nt, nt_guess
335 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
336 : LOGICAL :: all_conv, has_prev, stalled
337 2 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
338 : REAL(KIND=dp) :: dev, t_iter
339 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, diag, diag_sq, omega, omega_prev, &
340 2 : omega_sq, res, theta, theta_sub
341 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_restart, coef_restart_cand, &
342 2 : coef_ritz, coef_x, coef_x_omega, &
343 2 : coef_y, red_block
344 : TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
345 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
346 : TYPE(cp_fm_type) :: fm_metric, fm_metric_vec, fm_MW, &
347 : fm_red_M, fm_scratch, fm_V, fm_W, &
348 : fm_work
349 : TYPE(mp_para_env_type), POINTER :: para_env
350 :
351 2 : CALL timeset(routineN, handle)
352 2 : CALL cite_reference(Davidson1975)
353 2 : CALL cite_reference(Vecharynski2017)
354 2 : CALL cite_reference(Bai2012)
355 2 : CALL cite_reference(Stathopoulos1998)
356 2 : CALL cite_reference(Fukaya2014)
357 :
358 2 : abba_status = abba_ok
359 2 : ab_margin = 0.0_dp
360 2 : para_env => mv_env%para_env
361 : ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
362 2 : CPASSERT(mv_env%blacs_env%num_pe(2) == 1)
363 : ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
364 2 : NULLIFY (diag_blacs_env)
365 2 : CALL cp_blacs_env_create(diag_blacs_env, para_env)
366 2 : n_ov = mv_env%n_ov
367 2 : CALL davidson_sizes(mp2_env, n_ov, unit_nr, n_want, n_act, block_size)
368 :
369 6 : ALLOCATE (diag(n_ov))
370 2 : CALL bse_matvec_diagonal(mv_env, mp2_env%bse%preconditioner, diag)
371 :
372 2 : CALL subspace_ceiling(mp2_env, mv_env, driver_mk, n_act, block_size, unit_nr, m_max)
373 :
374 2 : CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
375 2 : CALL cp_fm_create(fm_V, fm_struct, name="fm_V_davidson", set_zero=.TRUE.)
376 2 : CALL cp_fm_create(fm_W, fm_struct, name="fm_W_davidson", set_zero=.TRUE.)
377 2 : CALL cp_fm_create(fm_MW, fm_struct, name="fm_MW_davidson", set_zero=.TRUE.)
378 2 : CALL cp_fm_struct_release(fm_struct)
379 2 : CALL bse_matvec_vector_struct(mv_env, MIN(2*n_act, m_max), fm_struct)
380 2 : CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.TRUE.)
381 2 : CALL cp_fm_struct_release(fm_struct)
382 2 : CALL bse_matvec_vector_struct(mv_env, block_size, fm_struct)
383 2 : CALL cp_fm_create(fm_scratch, fm_struct, name="fm_scratch_davidson", set_zero=.TRUE.)
384 2 : CALL cp_fm_struct_release(fm_struct)
385 2 : NULLIFY (fm_struct)
386 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
387 2 : nrow_global=m_max, ncol_global=m_max)
388 2 : CALL cp_fm_create(fm_red_M, fm_struct, name="fm_red_M_davidson", set_zero=.TRUE.)
389 2 : CALL cp_fm_struct_release(fm_struct)
390 2 : CALL cp_fm_get_info(fm_V, nrow_local=nrow_local)
391 :
392 : ! c_mk, the restart candidates and basis, then (Ω^k)^1/2 c_mk, (Ω^k)^-1/2 c_mk and Ω^k (Ω^k)^1/2 c_mk
393 20 : ALLOCATE (coef_ritz(m_max, n_act), coef_restart_cand(m_max, 2*n_act), coef_restart(m_max, 2*n_act))
394 14 : ALLOCATE (coef_x(m_max, n_act), coef_y(m_max, n_act), coef_x_omega(m_max, n_act))
395 : ALLOCATE (theta(m_max), omega(n_act), omega_prev(n_act), res(n_act), dE(n_act), conv(n_act), &
396 24 : need(n_act), selected(n_act))
397 2 : coef_restart_cand(:, :) = 0.0_dp
398 2 : omega_prev(:) = 0.0_dp
399 2 : has_prev = .FALSE.
400 2 : stalled = .FALSE.
401 2 : all_conv = .FALSE.
402 2 : n_restart = 0
403 2 : n_kernel = 0
404 2 : n_dependent = 0
405 :
406 : ! guess: the TDA vectors if given, then the eigenvectors of the exact guess block or unit
407 : ! vectors on the lowest entries of the diagonal
408 2 : n_seed = 0
409 2 : IF (PRESENT(fm_X_tda)) THEN
410 0 : CALL cp_fm_get_info(fm_X_tda, ncol_global=n_seed)
411 0 : n_seed = MIN(n_seed, n_act)
412 0 : CALL cp_fm_to_fm(fm_X_tda, fm_V, n_seed)
413 : END IF
414 2 : IF (mp2_env%bse%num_guess_transitions > 0) THEN
415 : CALL initial_guess_subblock(mv_env, mp2_env, diag, n_act, m_max, deg_thresh, fm_V, n_guess_cols, theta_sub, &
416 0 : unit_nr, n_seed, abba_status, ab_margin)
417 : ELSE
418 2 : CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_V, n_guess_cols, n_seed)
419 : END IF
420 : ! the corrections divide by d_ia^2 - (Ω^k)^2
421 4 : ALLOCATE (diag_sq(n_ov))
422 98 : diag_sq(:) = diag(:)**2
423 :
424 : ! the guess enters the basis in blocks of the scratch width
425 2 : m = 0
426 2 : guess_col = 1
427 4 : DO WHILE (guess_col <= n_guess_cols)
428 2 : nt_guess = MIN(block_size, n_guess_cols - guess_col + 1)
429 2 : nt = nt_guess
430 2 : IF (guess_col > m + 1) THEN
431 0 : CALL cp_fm_to_fm(fm_V, fm_V, nt, guess_col, m + 1)
432 : END IF
433 : CALL extend_basis_mk(mv_env, fm_V, fm_W, fm_MW, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
434 2 : ab_margin)
435 2 : IF (abba_status /= abba_ok) EXIT
436 2 : guess_col = guess_col + nt_guess
437 2 : m = m + nt
438 : END DO
439 :
440 2 : IF (abba_status == abba_ok) THEN
441 2 : fm_V%local_data(1:nrow_local, m + 1:n_guess_cols) = 0.0_dp
442 : ! Hr_mn = sum_ia W_ia,m (M W)_ia,n on the guess basis
443 2 : CALL reduced_gram_blocks(fm_W, m, fm_MW, 1, m, block_size, para_env, fm_red_M)
444 :
445 2 : CALL print_iteration_header('Block Davidson iterations for ABBA, (A+B)(A-B) x = E^2 x:', unit_nr)
446 :
447 16 : DO iter = 1, mp2_env%bse%max_iter
448 16 : t_iter = m_walltime()
449 :
450 : ! sum_n Hr_mn c_nk = θ_k c_mk with θ_k = (Ω^k)^2, ascending
451 16 : CALL solve_reduced(fm_red_M, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
452 :
453 : ! thick restart: V, W and MW are rotated alike by Q_pk in coef_restart, which keeps
454 : ! sum_ia V_ia,m W_ia,n = δ_mn, and Hr_kl <- sum_pq Q_pk Hr_pq Q_ql
455 16 : IF (m + block_size > m_max .AND. m > n_act .AND. m < n_ov) THEN
456 324 : coef_restart_cand(1:m, 1:n_act) = coef_ritz(1:m, 1:n_act)
457 2 : n_cand = n_act
458 2 : IF (has_prev) n_cand = 2*n_act
459 2 : CALL restart_basis(coef_restart_cand, m, n_cand, 2*n_act, coef_restart, k)
460 2 : CALL rotate_in_place(fm_V, m, coef_restart, k, fm_work, fm_W, fm_MW)
461 2 : CALL reduced_rotate(fm_red_M, m, coef_restart, k, para_env, diag_blacs_env)
462 2 : m = k
463 2 : n_restart = n_restart + 1
464 2 : CALL solve_reduced(fm_red_M, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
465 : END IF
466 :
467 : ! a Ritz value bounds an eigenvalue of (A+B)(A-B) from above
468 16 : IF (theta(1) <= 0.0_dp) THEN
469 : CALL cp_abort(__LOCATION__, &
470 : "BSE Davidson (MK_DAVIDSON): negative squared excitation energy. Matrix "// &
471 0 : "(A+B) is not positive definite, or the basis lost its orthonormality.")
472 : END IF
473 128 : omega(:) = SQRT(theta(1:n_act))
474 :
475 : ! r_ia,k = sum_m (M W)_ia,m c_mk (Ω^k)^-1/2 - Ω^k (Ω^k)^1/2 sum_m V_ia,m c_mk
476 128 : DO k = 1, n_act
477 1792 : coef_y(1:m, k) = coef_ritz(1:m, k)/SQRT(omega(k))
478 1808 : coef_x_omega(1:m, k) = coef_ritz(1:m, k)*omega(k)*SQRT(omega(k))
479 : END DO
480 1808 : CALL subspace_rotate(fm_MW, m, coef_y(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
481 1808 : CALL subspace_rotate(fm_V, m, coef_x_omega(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
482 16 : CALL column_norms(fm_work, 1, n_act, para_env, res)
483 :
484 : CALL convergence_step(mp2_env, omega, omega_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
485 16 : dE, conv, n_req, need)
486 38 : all_conv = .NOT. ANY(need)
487 16 : IF (all_conv) EXIT
488 :
489 14 : IF (iter == mp2_env%bse%max_iter) EXIT
490 :
491 : ! corrections t_ia,k = r_ia,k/(d_ia^2 - (Ω^k)^2) for the lowest unconverged roots
492 14 : CALL select_roots(mp2_env, need, res, MIN(block_size, m_max - m), selected, nt)
493 14 : CALL davidson_corrections(fm_work, selected(1:nt), theta, diag_sq, fm_V, m + 1)
494 : CALL extend_basis_mk(mv_env, fm_V, fm_W, fm_MW, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
495 14 : ab_margin)
496 14 : IF (abba_status /= abba_ok) EXIT
497 : ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
498 : ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
499 14 : IF (nt == 0 .AND. stalled) THEN
500 0 : CPABORT("BSE Davidson: correction vectors vanished before convergence")
501 : END IF
502 14 : stalled = nt == 0
503 :
504 14 : IF (nt > 0) THEN
505 : ! new columns Hr_mn = sum_ia W_ia,m (M W)_ia,n of the symmetric reduced matrix
506 56 : ALLOCATE (red_block(m + nt, nt))
507 14 : CALL subspace_gram(fm_W, 1, m + nt, fm_MW, m + 1, nt, para_env, red_block)
508 14 : CALL reduced_extend(fm_red_M, red_block, m, nt)
509 14 : DEALLOCATE (red_block)
510 : END IF
511 :
512 : ! the Ritz vectors of this iteration are the second block of restart candidates
513 2856 : coef_restart_cand(:, n_act + 1:2*n_act) = 0.0_dp
514 1540 : coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_ritz(1:m, 1:n_act)
515 112 : has_prev = .TRUE.
516 112 : omega_prev(:) = omega(:)
517 32 : m = m + nt
518 : END DO
519 : END IF
520 :
521 2 : IF (abba_status == abba_ok) THEN
522 2 : IF (.NOT. all_conv) CALL abort_unconverged(omega, res, conv, need, n_req, unit_nr)
523 2 : IF (mp2_env%bse%bse_debug_print) CALL print_tracked_roots(omega, res, n_want, n_act, unit_nr)
524 :
525 2 : IF (mp2_env%bse%bse_debug_print) THEN
526 0 : dev = orthonormality_deviation(fm_V, fm_W, m, block_size, para_env)
527 0 : IF (unit_nr > 0) THEN
528 0 : WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
529 0 : 'Max deviation of the basis from orthonormality', dev
530 : END IF
531 : END IF
532 :
533 : ! with sum_ia V_ia,m W_ia,n = δ_mn the smallest eigenvalue of A-B on the span of V is the
534 : ! inverse of the largest eigenvalue of the metric sum_ia V_ia,m V_ia,n
535 2 : NULLIFY (fm_struct)
536 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
537 2 : nrow_global=m, ncol_global=m)
538 2 : CALL cp_fm_create(fm_metric, fm_struct, name="bse_basis_metric")
539 2 : CALL cp_fm_create(fm_metric_vec, fm_struct, name="bse_basis_metric_vectors")
540 2 : CALL cp_fm_struct_release(fm_struct)
541 2 : CALL reduced_gram_blocks(fm_V, m, fm_V, 1, m, block_size, para_env, fm_metric)
542 6 : ALLOCATE (omega_sq(m))
543 2 : omega_sq(:) = 0.0_dp
544 2 : CALL choose_eigv_solver(fm_metric, fm_metric_vec, omega_sq)
545 2 : ab_margin = 1.0_dp/omega_sq(m)
546 2 : DEALLOCATE (omega_sq)
547 2 : CALL cp_fm_release(fm_metric)
548 2 : CALL cp_fm_release(fm_metric_vec)
549 :
550 2 : IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, omega, n_req, unit_nr)
551 2 : CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, n_dependent)
552 :
553 : ! x_ia^n = (Ω^n)^1/2 sum_m V_ia,m c_mn, y_ia^n = (Ω^n)^-1/2 sum_m W_ia,m c_mn,
554 : ! X_ia^n = (y_ia^n + x_ia^n)/2, Y_ia^n = (y_ia^n - x_ia^n)/2
555 6 : ALLOCATE (exc_ens(n_want))
556 8 : exc_ens(:) = omega(1:n_want)
557 8 : DO k = 1, n_want
558 114 : coef_x(1:m, k) = coef_ritz(1:m, k)*SQRT(omega(k))
559 116 : coef_y(1:m, k) = coef_ritz(1:m, k)/SQRT(omega(k))
560 : END DO
561 2 : CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
562 2 : CALL cp_fm_create(fm_X, fm_struct, name="fm_X_davidson", set_zero=.TRUE.)
563 2 : CALL cp_fm_create(fm_Y, fm_struct, name="fm_Y_davidson", set_zero=.TRUE.)
564 2 : CALL cp_fm_struct_release(fm_struct)
565 116 : CALL subspace_rotate(fm_W, m, coef_y(1:m, 1:n_want), n_want, fm_X, 0.5_dp, 0.0_dp)
566 116 : CALL subspace_rotate(fm_V, m, coef_x(1:m, 1:n_want), n_want, fm_X, 0.5_dp, 1.0_dp)
567 116 : CALL subspace_rotate(fm_W, m, coef_y(1:m, 1:n_want), n_want, fm_Y, 0.5_dp, 0.0_dp)
568 118 : CALL subspace_rotate(fm_V, m, coef_x(1:m, 1:n_want), n_want, fm_Y, -0.5_dp, 1.0_dp)
569 : END IF
570 :
571 2 : CALL cp_fm_release(fm_V)
572 2 : CALL cp_fm_release(fm_W)
573 2 : CALL cp_fm_release(fm_MW)
574 2 : CALL cp_fm_release(fm_work)
575 2 : CALL cp_fm_release(fm_scratch)
576 2 : CALL cp_fm_release(fm_red_M)
577 2 : CALL cp_blacs_env_release(diag_blacs_env)
578 0 : DEALLOCATE (coef_ritz, coef_restart_cand, coef_x, coef_y, coef_x_omega, coef_restart, theta, omega, omega_prev, &
579 2 : res, dE, conv, need, selected, diag, diag_sq)
580 2 : IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
581 :
582 2 : CALL timestop(handle)
583 :
584 8 : END SUBROUTINE bse_davidson_abba_mk
585 :
586 : ! **************************************************************************************************
587 : !> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem
588 : !> in the paired form of Olsen and Stratmann: with x^n = X^n - Y^n, y^n = X^n + Y^n, M = A+B and
589 : !> K = A-B, sum_jb M_ia,jb y_jb^n = Ω^n x_ia^n and sum_jb K_ia,jb x_jb^n = Ω^n y_ia^n. One
590 : !> orthonormal basis b_ia,m carries both vectors: Mr_mn = sum_ia b_ia,m (M b)_ia,n and
591 : !> Kr_mn = sum_ia b_ia,m (K b)_ia,n are solved by solve_reduced_paired for Ω^k, R_mk and L_mk with
592 : !> y_ia^k = sum_m b_ia,m R_mk, x_ia^k = sum_m b_ia,m L_mk and sum_ia x_ia^k y_ia^k = 1. The basis
593 : !> grows by the corrections t_ia,k = r_ia,k/(d_ia - Ω^k) of both residuals,
594 : !> rR_ia,k = sum_m (M b)_ia,m R_mk - Ω^k x_ia^k and rL_ia,k = sum_m (K b)_ia,m L_mk - Ω^k y_ia^k,
595 : !> d_ia chosen by PRECONDITIONER; the rR block enters before the rL block, which is nearly
596 : !> parallel to it for a small B. Arrays: Mr_mn and Kr_mn in fm_red_M and fm_red_K on the process
597 : !> grid, R_mk in coef_right, L_mk in coef_left
598 : !> \param mv_env slabs, prefactors and transition energies of the matrix-free A and B
599 : !> \param mp2_env mp2_env%bse carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
600 : !> \param unit_nr output unit, positive on the writing rank only
601 : !> \param exc_ens Ω^n in Hartree
602 : !> \param fm_X X_ia^n on the process grid of mv_env, column n; created only for abba_ok
603 : !> \param fm_Y Y_ia^n, as fm_X
604 : !> \param abba_status abba_indefinite if A-B is not positive definite on the trial space
605 : !> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
606 : !> \param fm_X_tda converged TDA vectors as initial guess, x = y = X
607 : ! **************************************************************************************************
608 2 : SUBROUTINE bse_davidson_abba_os(mv_env, mp2_env, unit_nr, exc_ens, fm_X, fm_Y, abba_status, &
609 : ab_margin, fm_X_tda)
610 :
611 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
612 : TYPE(mp2_type), INTENT(IN) :: mp2_env
613 : INTEGER, INTENT(IN) :: unit_nr
614 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
615 : INTENT(OUT) :: exc_ens
616 : TYPE(cp_fm_type), INTENT(OUT) :: fm_X, fm_Y
617 : INTEGER, INTENT(OUT) :: abba_status
618 : REAL(KIND=dp), INTENT(OUT) :: ab_margin
619 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_X_tda
620 :
621 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_davidson_abba_os'
622 :
623 : INTEGER :: block_size, guess_col, handle, iter, k, m, m_max, n_act, n_cand, n_dependent, &
624 : n_guess_cols, n_kernel, n_ov, n_pair, n_req, n_restart, n_seed, n_want, nrow_local, nt, &
625 : nt_guess, nt_L, nt_R, nt_sel
626 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
627 : LOGICAL :: all_conv, has_prev, stalled
628 2 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
629 : REAL(KIND=dp) :: dev, kr_min, t_iter
630 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dE, diag, omega, omega_prev, res, res_L, &
631 2 : theta_sub
632 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_left, coef_restart, &
633 : coef_restart_cand, coef_right, &
634 2 : coef_work, red_block_K, red_block_M
635 : TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
636 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
637 : TYPE(cp_fm_type) :: fm_b, fm_Kb, fm_Mb, fm_red_K, fm_red_M, &
638 : fm_scratch, fm_work
639 : TYPE(mp_para_env_type), POINTER :: para_env
640 :
641 2 : CALL timeset(routineN, handle)
642 2 : CALL cite_reference(Davidson1975)
643 2 : CALL cite_reference(Olsen1988)
644 2 : CALL cite_reference(Stratmann1998)
645 2 : CALL cite_reference(Bai2012)
646 2 : CALL cite_reference(Stathopoulos1998)
647 2 : CALL cite_reference(Fukaya2014)
648 :
649 2 : abba_status = abba_ok
650 2 : ab_margin = 0.0_dp
651 2 : para_env => mv_env%para_env
652 : ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
653 2 : CPASSERT(mv_env%blacs_env%num_pe(2) == 1)
654 : ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
655 2 : NULLIFY (diag_blacs_env)
656 2 : CALL cp_blacs_env_create(diag_blacs_env, para_env)
657 2 : n_ov = mv_env%n_ov
658 2 : CALL davidson_sizes(mp2_env, n_ov, unit_nr, n_want, n_act, block_size)
659 :
660 6 : ALLOCATE (diag(n_ov))
661 2 : CALL bse_matvec_diagonal(mv_env, mp2_env%bse%preconditioner, diag)
662 :
663 2 : CALL subspace_ceiling(mp2_env, mv_env, driver_os, n_act, block_size, unit_nr, m_max)
664 :
665 2 : CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
666 2 : CALL cp_fm_create(fm_b, fm_struct, name="fm_b_davidson", set_zero=.TRUE.)
667 2 : CALL cp_fm_create(fm_Mb, fm_struct, name="fm_Mb_davidson", set_zero=.TRUE.)
668 2 : CALL cp_fm_create(fm_Kb, fm_struct, name="fm_Kb_davidson", set_zero=.TRUE.)
669 2 : CALL cp_fm_struct_release(fm_struct)
670 : ! both residual blocks, or the restart basis
671 2 : CALL bse_matvec_vector_struct(mv_env, MAX(2*n_act, MIN(4*n_act, m_max)), fm_struct)
672 2 : CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.TRUE.)
673 2 : CALL cp_fm_struct_release(fm_struct)
674 2 : CALL bse_matvec_vector_struct(mv_env, 2*block_size, fm_struct)
675 2 : CALL cp_fm_create(fm_scratch, fm_struct, name="fm_scratch_davidson", set_zero=.TRUE.)
676 2 : CALL cp_fm_struct_release(fm_struct)
677 2 : NULLIFY (fm_struct)
678 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
679 2 : nrow_global=m_max, ncol_global=m_max)
680 2 : CALL cp_fm_create(fm_red_M, fm_struct, name="fm_red_M_davidson", set_zero=.TRUE.)
681 2 : CALL cp_fm_create(fm_red_K, fm_struct, name="fm_red_K_davidson", set_zero=.TRUE.)
682 2 : CALL cp_fm_struct_release(fm_struct)
683 2 : CALL cp_fm_get_info(fm_b, nrow_local=nrow_local)
684 :
685 : ! R_mk, L_mk and a work block, then the restart candidates (current and previous pairs) and the basis Q_pk
686 16 : ALLOCATE (coef_right(m_max, n_act), coef_left(m_max, n_act), coef_work(m_max, n_act))
687 14 : ALLOCATE (coef_restart_cand(m_max, 4*n_act), coef_restart(m_max, 4*n_act))
688 : ALLOCATE (omega(n_act), omega_prev(n_act), res(n_act), res_L(n_act), dE(n_act), conv(n_act), &
689 22 : need(n_act), selected(n_act))
690 2 : coef_restart_cand(:, :) = 0.0_dp
691 2 : omega_prev(:) = 0.0_dp
692 2 : kr_min = 0.0_dp
693 2 : has_prev = .FALSE.
694 2 : stalled = .FALSE.
695 2 : all_conv = .FALSE.
696 2 : n_restart = 0
697 2 : n_kernel = 0
698 2 : n_dependent = 0
699 :
700 : ! guess: the TDA vectors if given (each a pair by itself), then the pairs x_k, y_k of the exact
701 : ! guess block or unit vectors on the lowest entries of the diagonal
702 2 : n_seed = 0
703 2 : IF (PRESENT(fm_X_tda)) THEN
704 0 : CALL cp_fm_get_info(fm_X_tda, ncol_global=n_seed)
705 0 : n_seed = MIN(n_seed, n_act)
706 0 : CALL cp_fm_to_fm(fm_X_tda, fm_b, n_seed)
707 : END IF
708 2 : n_pair = 0
709 2 : IF (mp2_env%bse%num_guess_transitions > 0) THEN
710 : CALL initial_guess_subblock(mv_env, mp2_env, diag, n_act, m_max, deg_thresh, fm_b, n_guess_cols, theta_sub, &
711 0 : unit_nr, n_seed, abba_status, ab_margin, n_pair)
712 : ELSE
713 2 : CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_b, n_guess_cols, n_seed)
714 : END IF
715 2 : n_guess_cols = n_guess_cols + n_pair
716 :
717 : ! the guess enters the basis in blocks of the scratch width
718 2 : m = 0
719 2 : guess_col = 1
720 4 : DO WHILE (guess_col <= n_guess_cols)
721 2 : nt_guess = MIN(2*block_size, n_guess_cols - guess_col + 1)
722 2 : nt = nt_guess
723 2 : IF (guess_col > m + 1) THEN
724 0 : CALL cp_fm_to_fm(fm_b, fm_b, nt, guess_col, m + 1)
725 : END IF
726 2 : CALL extend_basis_os(mv_env, fm_b, fm_Mb, fm_Kb, fm_scratch, m, nt, n_kernel, n_dependent)
727 2 : guess_col = guess_col + nt_guess
728 2 : m = m + nt
729 : END DO
730 :
731 2 : IF (abba_status == abba_ok) THEN
732 2 : fm_b%local_data(1:nrow_local, m + 1:n_guess_cols) = 0.0_dp
733 : ! Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n on the guess basis
734 2 : CALL reduced_gram_blocks(fm_b, m, fm_Mb, 1, m, block_size, para_env, fm_red_M)
735 2 : CALL reduced_gram_blocks(fm_b, m, fm_Kb, 1, m, block_size, para_env, fm_red_K)
736 :
737 : CALL print_iteration_header('Block Davidson iterations for ABBA, Olsen-Stratmann paired subspace:', &
738 2 : unit_nr)
739 :
740 18 : DO iter = 1, mp2_env%bse%max_iter
741 18 : t_iter = m_walltime()
742 :
743 : ! Ω^k ascending with the coefficients R_mk of y^k and L_mk of x^k, sum_m R_mk L_mk = 1
744 : CALL solve_reduced_paired(fm_red_M, fm_red_K, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
745 18 : abba_status, diag_blacs_env)
746 : ! the smallest eigenvalue of A-B on the span of b is that of the last Kr
747 18 : ab_margin = kr_min
748 18 : IF (abba_status /= abba_ok) EXIT
749 :
750 : ! thick restart: b, Mb and Kb are rotated alike; the current pairs come first, so that
751 : ! they survive the cap that leaves room for one block of corrections; with Q_pk in
752 : ! coef_restart, b_ia,k <- sum_p b_ia,p Q_pk and Mr_kl <- sum_pq Q_pk Mr_pq Q_ql, Kr alike
753 18 : IF (m + 2*block_size > m_max .AND. m > 2*n_act .AND. m < n_ov) THEN
754 2214 : coef_restart_cand(1:m, 1:n_act) = coef_right(1:m, :)
755 2214 : coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_left(1:m, :)
756 16 : n_cand = 2*n_act
757 16 : IF (has_prev) n_cand = 4*n_act
758 16 : CALL restart_basis(coef_restart_cand, m, n_cand, m_max - 2*block_size, coef_restart, k)
759 16 : CALL rotate_in_place(fm_b, m, coef_restart, k, fm_work, fm_Mb, fm_Kb)
760 16 : CALL reduced_rotate(fm_red_M, m, coef_restart, k, para_env, diag_blacs_env)
761 16 : CALL reduced_rotate(fm_red_K, m, coef_restart, k, para_env, diag_blacs_env)
762 16 : m = k
763 16 : n_restart = n_restart + 1
764 : CALL solve_reduced_paired(fm_red_M, fm_red_K, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
765 16 : abba_status, diag_blacs_env)
766 16 : ab_margin = kr_min
767 16 : IF (abba_status /= abba_ok) EXIT
768 : END IF
769 :
770 : ! rR_ia,k = sum_m (M b)_ia,m R_mk - Ω^k sum_m b_ia,m L_mk into the columns 1..n_act and
771 : ! rL_ia,k = sum_m (K b)_ia,m L_mk - Ω^k sum_m b_ia,m R_mk into the columns n_act+1..2 n_act
772 144 : DO k = 1, n_act
773 1810 : coef_work(1:m, k) = coef_left(1:m, k)*omega(k)
774 : END DO
775 1810 : CALL subspace_rotate(fm_Mb, m, coef_right(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
776 1810 : CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
777 144 : DO k = 1, n_act
778 1810 : coef_work(1:m, k) = coef_right(1:m, k)*omega(k)
779 : END DO
780 1810 : CALL subspace_rotate(fm_Kb, m, coef_left(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp, n_act + 1)
781 1810 : CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp, n_act + 1)
782 18 : CALL column_norms(fm_work, 1, n_act, para_env, res)
783 18 : CALL column_norms(fm_work, n_act + 1, n_act, para_env, res_L)
784 144 : res(:) = MAX(res(:), res_L(:))
785 :
786 : CALL convergence_step(mp2_env, omega, omega_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
787 18 : dE, conv, n_req, need)
788 32 : all_conv = .NOT. ANY(need)
789 18 : IF (all_conv) EXIT
790 :
791 16 : IF (iter == mp2_env%bse%max_iter) EXIT
792 :
793 : ! corrections t_ia,k = r_ia,k/(d_ia - Ω^k) of the lowest unconverged roots: the rR block
794 : ! against the basis, then the rL block of the same roots against the extended basis;
795 : ! both fit the room, which is odd only when the ceiling is the number of transitions
796 16 : CALL select_roots(mp2_env, need, res, MIN(block_size, (m_max - m + 1)/2), selected, nt_sel)
797 16 : nt_R = nt_sel
798 16 : CALL davidson_corrections(fm_work, selected(1:nt_R), omega, diag, fm_b, m + 1)
799 16 : CALL extend_basis_os(mv_env, fm_b, fm_Mb, fm_Kb, fm_scratch, m, nt_R, n_kernel, n_dependent)
800 16 : nt_L = MIN(nt_sel, m_max - m - nt_R)
801 16 : CALL davidson_corrections(fm_work, selected(1:nt_L), omega, diag, fm_b, m + nt_R + 1, n_act)
802 16 : CALL extend_basis_os(mv_env, fm_b, fm_Mb, fm_Kb, fm_scratch, m + nt_R, nt_L, n_kernel, n_dependent)
803 16 : nt = nt_R + nt_L
804 : ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
805 : ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
806 16 : IF (nt == 0 .AND. stalled) THEN
807 0 : CPABORT("BSE Davidson: correction vectors vanished before convergence")
808 : END IF
809 16 : stalled = nt == 0
810 :
811 16 : IF (nt > 0) THEN
812 : ! new columns n of Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n
813 96 : ALLOCATE (red_block_M(m + nt, nt), red_block_K(m + nt, nt))
814 16 : CALL subspace_gram(fm_b, 1, m + nt, fm_Mb, m + 1, nt, para_env, red_block_M)
815 16 : CALL subspace_gram(fm_b, 1, m + nt, fm_Kb, m + 1, nt, para_env, red_block_K)
816 16 : CALL reduced_extend(fm_red_M, red_block_M, m, nt)
817 16 : CALL reduced_extend(fm_red_K, red_block_K, m, nt)
818 16 : DEALLOCATE (red_block_M, red_block_K)
819 : END IF
820 :
821 : ! the pairs of this iteration are the second block of restart candidates
822 6512 : coef_restart_cand(:, 2*n_act + 1:4*n_act) = 0.0_dp
823 1598 : coef_restart_cand(1:m, 2*n_act + 1:3*n_act) = coef_right(1:m, :)
824 1598 : coef_restart_cand(1:m, 3*n_act + 1:4*n_act) = coef_left(1:m, :)
825 128 : has_prev = .TRUE.
826 128 : omega_prev(:) = omega(:)
827 54 : m = m + nt
828 : END DO
829 : END IF
830 :
831 2 : IF (abba_status == abba_ok) THEN
832 2 : IF (.NOT. all_conv) CALL abort_unconverged(omega, res, conv, need, n_req, unit_nr)
833 2 : IF (mp2_env%bse%bse_debug_print) CALL print_tracked_roots(omega, res, n_want, n_act, unit_nr)
834 :
835 2 : IF (mp2_env%bse%bse_debug_print) THEN
836 0 : dev = orthonormality_deviation(fm_b, fm_b, m, block_size, para_env)
837 0 : IF (unit_nr > 0) THEN
838 0 : WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
839 0 : 'Max deviation of the basis from orthonormality', dev
840 : END IF
841 : END IF
842 :
843 2 : IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, omega, n_req, unit_nr)
844 2 : CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, n_dependent)
845 :
846 : ! X_ia^n = sum_m b_ia,m (R_mn + L_mn)/2, Y_ia^n = sum_m b_ia,m (R_mn - L_mn)/2
847 6 : ALLOCATE (exc_ens(n_want))
848 8 : exc_ens(:) = omega(1:n_want)
849 2 : CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
850 2 : CALL cp_fm_create(fm_X, fm_struct, name="fm_X_davidson", set_zero=.TRUE.)
851 2 : CALL cp_fm_create(fm_Y, fm_struct, name="fm_Y_davidson", set_zero=.TRUE.)
852 2 : CALL cp_fm_struct_release(fm_struct)
853 92 : coef_work(1:m, 1:n_want) = 0.5_dp*(coef_right(1:m, 1:n_want) + coef_left(1:m, 1:n_want))
854 92 : CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_want), n_want, fm_X, 1.0_dp, 0.0_dp)
855 92 : coef_work(1:m, 1:n_want) = 0.5_dp*(coef_right(1:m, 1:n_want) - coef_left(1:m, 1:n_want))
856 92 : CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_want), n_want, fm_Y, 1.0_dp, 0.0_dp)
857 : END IF
858 :
859 2 : CALL cp_fm_release(fm_b)
860 2 : CALL cp_fm_release(fm_Mb)
861 2 : CALL cp_fm_release(fm_Kb)
862 2 : CALL cp_fm_release(fm_work)
863 2 : CALL cp_fm_release(fm_scratch)
864 2 : CALL cp_fm_release(fm_red_M)
865 2 : CALL cp_fm_release(fm_red_K)
866 2 : CALL cp_blacs_env_release(diag_blacs_env)
867 0 : DEALLOCATE (coef_right, coef_left, coef_work, coef_restart_cand, coef_restart, omega, omega_prev, res, res_L, &
868 2 : dE, conv, need, selected, diag)
869 2 : IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
870 :
871 2 : CALL timestop(handle)
872 :
873 8 : END SUBROUTINE bse_davidson_abba_os
874 :
875 : ! **************************************************************************************************
876 : !> \brief Makes the nt columns T that follow column m of V orthonormal in the (A-B) inner product and
877 : !> stores their products: T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb W_jb,m T_jb,k, Euclidean
878 : !> orthonormalisation of T, W_ia,k = sum_jb (A-B)_ia,jb T_jb,k, T_ia,k <- sum_l T_ia,l (R^-1)_lk and
879 : !> W alike with sum_p R_pk R_pl = sum_ia T_ia,k W_ia,l, and (M W)_ia,k = sum_jb (A+B)_ia,jb W_jb,k
880 : !> \param mv_env the matrix-free A and B
881 : !> \param fm_V basis; columns 1..m are kept, m+1..m+nt hold T on entry and the kept columns on exit
882 : !> \param fm_W (A-B) V, extended in the same columns
883 : !> \param fm_MW (A+B) W, extended in the same columns
884 : !> \param fm_scratch scratch of at least nt columns
885 : !> \param m columns of the basis before the extension
886 : !> \param nt on exit the number of columns kept
887 : !> \param n_kernel counter of the kernel applications
888 : !> \param n_dependent counter of the correction vectors dropped as dependent
889 : !> \param abba_status abba_indefinite if the Gram matrix of T in the (A-B) inner product is not positive definite
890 : !> \param ab_margin smallest eigenvalue of T^T (A-B) T if that matrix is not positive definite
891 : ! **************************************************************************************************
892 16 : SUBROUTINE extend_basis_mk(mv_env, fm_V, fm_W, fm_MW, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
893 : ab_margin)
894 :
895 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
896 : TYPE(cp_fm_type), INTENT(IN) :: fm_V, fm_W, fm_MW, fm_scratch
897 : INTEGER, INTENT(IN) :: m
898 : INTEGER, INTENT(INOUT) :: nt, n_kernel, n_dependent, abba_status
899 : REAL(KIND=dp), INTENT(INOUT) :: ab_margin
900 :
901 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extend_basis_mk'
902 :
903 : INTEGER :: handle
904 :
905 16 : CALL timeset(routineN, handle)
906 :
907 : ! T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb W_jb,m T_jb,k, then Euclidean orthonormalisation
908 16 : CALL extend_orthonormal(fm_V, m, nt, mv_env%para_env, n_dependent, fm_W)
909 :
910 16 : IF (nt > 0) THEN
911 : ! W_ia,k = sum_jb (A-B)_ia,jb T_jb,k: A T into W, B T into the scratch
912 16 : CALL bse_matvec_apply(mv_env, fm_V, m + 1, nt, fm_W, fm_scratch, 1)
913 16 : CALL columns_axpy(-1.0_dp, fm_scratch, 1, fm_W, m + 1, nt)
914 16 : n_kernel = n_kernel + nt
915 :
916 : ! T_ia,k <- sum_l T_ia,l (R^-1)_lk and W alike, with sum_p R_pk R_pl = sum_ia T_ia,k W_ia,l
917 16 : CALL cholqr2_metric(fm_V, fm_W, m + 1, nt, mv_env%para_env, abba_status, ab_margin)
918 :
919 16 : IF (abba_status == abba_ok) THEN
920 : ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
921 16 : IF (extension_deviation(fm_V, fm_W, m, nt, mv_env%para_env) > ortho_tol) THEN
922 0 : CPABORT("BSE Davidson: the basis lost orthonormality in the (A-B) inner product")
923 : END IF
924 : ! (M W)_ia,k = sum_jb (A+B)_ia,jb W_jb,k: A W into MW, B W into the scratch
925 16 : CALL bse_matvec_apply(mv_env, fm_W, m + 1, nt, fm_MW, fm_scratch, 1)
926 16 : CALL columns_axpy(1.0_dp, fm_scratch, 1, fm_MW, m + 1, nt)
927 16 : n_kernel = n_kernel + nt
928 : END IF
929 : END IF
930 :
931 16 : CALL timestop(handle)
932 :
933 16 : END SUBROUTINE extend_basis_mk
934 :
935 : ! **************************************************************************************************
936 : !> \brief Two Cholesky QR passes in the inner product of a symmetric K with the products KT at hand:
937 : !> G_kl = sum_ia T_ia,k (KT)_ia,l = sum_p R_pk R_pl, T_ia,k <- sum_l T_ia,l (R^-1)_lk, KT alike. A failed
938 : !> factorization of a G with a negative eigenvalue reports K as indefinite, any other failure
939 : !> is repeated once with a diagonal shift and a third pass. The Gram is the only communication:
940 : !> R is factorised on one rank and broadcast, the triangular solve runs on the local rows
941 : !> \param fm_T the block to orthonormalise, columns first_col .. first_col+nt-1
942 : !> \param fm_KT K T in the same columns, transformed alike
943 : !> \param first_col first column of the block
944 : !> \param nt number of columns in the block
945 : !> \param para_env communicator of the row-distributed matrices
946 : !> \param abba_status abba_indefinite if gram has a negative eigenvalue, else untouched
947 : !> \param eig_min smallest eigenvalue of gram, set if it is negative
948 : ! **************************************************************************************************
949 16 : SUBROUTINE cholqr2_metric(fm_T, fm_KT, first_col, nt, para_env, abba_status, eig_min)
950 :
951 : TYPE(cp_fm_type), INTENT(IN) :: fm_T, fm_KT
952 : INTEGER, INTENT(IN) :: first_col, nt
953 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
954 : INTEGER, INTENT(INOUT) :: abba_status
955 : REAL(KIND=dp), INTENT(INOUT) :: eig_min
956 :
957 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholqr2_metric'
958 :
959 : INTEGER :: handle, info, ipass, k, n_ov, n_pass, &
960 : nrow_local
961 : REAL(KIND=dp) :: shift
962 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
963 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chol_factor, gram, gram_raw
964 :
965 16 : CALL timeset(routineN, handle)
966 16 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local, nrow_global=n_ov)
967 160 : ALLOCATE (gram_raw(nt, nt), gram(nt, nt), chol_factor(nt, nt), eigval(nt))
968 :
969 16 : n_pass = 2
970 16 : ipass = 0
971 48 : DO WHILE (ipass < n_pass)
972 32 : ipass = ipass + 1
973 : ! G_kl = (sum_ia T_ia,k (KT)_ia,l + sum_ia T_ia,l (KT)_ia,k)/2 = sum_p R_pk R_pl
974 32 : CALL subspace_gram(fm_T, first_col, nt, fm_KT, first_col, nt, para_env, gram_raw)
975 608 : gram(:, :) = 0.5_dp*(gram_raw(:, :) + TRANSPOSE(gram_raw(:, :)))
976 608 : chol_factor(:, :) = gram(:, :)
977 32 : info = 0
978 32 : IF (para_env%is_source()) CALL DPOTRF('U', nt, chol_factor, nt, info)
979 32 : CALL para_env%bcast(info)
980 32 : IF (info /= 0) THEN
981 0 : IF (n_pass == 3) THEN
982 : CALL cp_abort(__LOCATION__, &
983 0 : "BSE Davidson: orthonormalisation in the (A-B) inner product broke down")
984 : END IF
985 0 : chol_factor(:, :) = gram(:, :)
986 0 : eigval(:) = 0.0_dp
987 0 : IF (para_env%is_source()) CALL diamat_all(chol_factor, eigval)
988 0 : CALL para_env%bcast(eigval)
989 0 : IF (eigval(1) < 0.0_dp) THEN
990 0 : abba_status = abba_indefinite
991 0 : eig_min = eigval(1)
992 0 : EXIT
993 : END IF
994 : ! a positive G that failed numerically: factorise G_kl + shift δ_kl and pass once more, the
995 : ! shift of shifted Cholesky QR (Fukaya et al. 2020), 11 (N_ov nt + nt (nt+1)) EPSILON |G|_F
996 0 : n_pass = 3
997 : shift = 11.0_dp*(REAL(n_ov, dp)*REAL(nt, dp) + REAL(nt, dp)*REAL(nt + 1, dp))* &
998 0 : EPSILON(1.0_dp)*NORM2(gram)
999 0 : chol_factor(:, :) = gram(:, :)
1000 0 : DO k = 1, nt
1001 0 : chol_factor(k, k) = chol_factor(k, k) + shift
1002 : END DO
1003 0 : IF (para_env%is_source()) CALL DPOTRF('U', nt, chol_factor, nt, info)
1004 0 : CALL para_env%bcast(info)
1005 0 : IF (info /= 0) THEN
1006 : CALL cp_abort(__LOCATION__, &
1007 0 : "BSE Davidson: orthonormalisation in the (A-B) inner product broke down")
1008 : END IF
1009 : END IF
1010 : ! T_ia,k <- sum_l T_ia,l (R^-1)_lk and KT alike, on the local rows
1011 32 : CALL para_env%bcast(chol_factor)
1012 48 : IF (nrow_local > 0) THEN
1013 : CALL DTRSM('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
1014 32 : fm_T%local_data(:, first_col:first_col + nt - 1), SIZE(fm_T%local_data, 1))
1015 : CALL DTRSM('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
1016 32 : fm_KT%local_data(:, first_col:first_col + nt - 1), SIZE(fm_KT%local_data, 1))
1017 : END IF
1018 : END DO
1019 :
1020 16 : DEALLOCATE (gram_raw, gram, chol_factor, eigval)
1021 16 : CALL timestop(handle)
1022 :
1023 16 : END SUBROUTINE cholqr2_metric
1024 :
1025 : ! **************************************************************************************************
1026 : !> \brief Paired reduced problem on the orthonormal basis b_ia,m,
1027 : !> Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n,
1028 : !> sum_n Mr_mn R_nk = Ω_k L_mk and sum_n Kr_mn L_nk = Ω_k R_mk,
1029 : !> solved as in initial_guess_subblock through the symmetric form:
1030 : !> Kr_mn = sum_p U_mp κ_p U_np and S_mn = sum_p U_mp √κ_p U_np, the square root of Kr,
1031 : !> sum_nop S_mn Mr_no S_op T_pk = λ_k T_mk with Ω_k = √λ_k,
1032 : !> R_mk = sum_n S_mn T_nk/√Ω_k and L_mk = sum_n Mr_mn R_nk/Ω_k.
1033 : !> κ_1 <= 0 reports K as indefinite on the span of b.
1034 : !> Every m x m step runs on the grid of blacs_env: the two eigensolves through
1035 : !> choose_eigv_solver and the four products through parallel_gemm; only the pair
1036 : !> coefficients come back replicated, R_mk in coef_right and L_mk in coef_left
1037 : !> \param fm_red_M Mr in its live block, on the grid of blacs_env
1038 : !> \param fm_red_K Kr in its live block, on the grid of blacs_env
1039 : !> \param m dimension of the live block
1040 : !> \param n_act number of pairs formed
1041 : !> \param para_env communicator of the grid
1042 : !> \param omega Ω_k in Hartree
1043 : !> \param coef_right coefficients of y_k in column k, zero below row m
1044 : !> \param coef_left coefficients of x_k in column k, zero below row m
1045 : !> \param kr_min κ_1, the smallest eigenvalue of Kr
1046 : !> \param abba_status abba_indefinite when κ_1 <= 0; omega, coef_right and coef_left are then untouched
1047 : !> \param blacs_env grid the reduced problem is solved on
1048 : ! **************************************************************************************************
1049 34 : SUBROUTINE solve_reduced_paired(fm_red_M, fm_red_K, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
1050 : abba_status, blacs_env)
1051 :
1052 : TYPE(cp_fm_type), INTENT(IN) :: fm_red_M, fm_red_K
1053 : INTEGER, INTENT(IN) :: m, n_act
1054 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
1055 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: omega
1056 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: coef_right, coef_left
1057 : REAL(KIND=dp), INTENT(OUT) :: kr_min
1058 : INTEGER, INTENT(INOUT) :: abba_status
1059 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
1060 :
1061 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_reduced_paired'
1062 :
1063 : INTEGER :: handle
1064 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: kappa, lambda
1065 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1066 : TYPE(cp_fm_type) :: fm_H, fm_H_eigvec, fm_K_eigvec, &
1067 : fm_K_sqrt, fm_L, fm_Mr, fm_R, fm_work
1068 :
1069 34 : CALL timeset(routineN, handle)
1070 : ! n_act pairs are extracted from an m-dimensional basis
1071 34 : CPASSERT(n_act <= m)
1072 :
1073 : ! Kr_mn = sum_p U_mp κ_p U_np, S_mn = sum_p U_mp √κ_p U_np, H_mp = sum_no S_mn Mr_no S_op, T its eigenvectors:
1074 : ! κ, U in fm_K_eigvec, S in fm_K_sqrt, the symmetrised Mr, H, T in fm_H_eigvec and a work matrix,
1075 : ! all m x m on the grid of blacs_env
1076 136 : ALLOCATE (kappa(m), lambda(m))
1077 34 : NULLIFY (fm_struct)
1078 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1079 34 : nrow_global=m, ncol_global=m)
1080 34 : CALL cp_fm_create(fm_K_eigvec, fm_struct, name="bse_paired_U")
1081 34 : CALL cp_fm_create(fm_K_sqrt, fm_struct, name="bse_paired_S")
1082 34 : CALL cp_fm_create(fm_Mr, fm_struct, name="bse_paired_Mr")
1083 34 : CALL cp_fm_create(fm_H, fm_struct, name="bse_paired_H")
1084 34 : CALL cp_fm_create(fm_H_eigvec, fm_struct, name="bse_paired_T")
1085 34 : CALL cp_fm_create(fm_work, fm_struct, name="bse_paired_work")
1086 34 : CALL cp_fm_struct_release(fm_struct)
1087 :
1088 : ! Kr_mn = sum_p U_mp κ_p U_np on the symmetrised live block, (Kr_mn + Kr_nm)/2
1089 34 : CALL cp_fm_to_fm_submat(fm_red_K, fm_work, m, m, 1, 1, 1, 1)
1090 34 : CALL symmetrise_in_place(fm_work, fm_K_eigvec)
1091 34 : kappa(:) = 0.0_dp
1092 34 : CALL choose_eigv_solver(fm_work, fm_K_eigvec, kappa)
1093 34 : kr_min = kappa(1)
1094 :
1095 34 : IF (kappa(1) <= 0.0_dp) THEN
1096 0 : abba_status = abba_indefinite
1097 : ELSE
1098 : ! S_mn = sum_p U_mp √κ_p U_np
1099 34 : CALL cp_fm_to_fm(fm_K_eigvec, fm_work)
1100 570 : CALL cp_fm_column_scale(fm_work, SQRT(kappa))
1101 34 : CALL parallel_gemm('N', 'T', m, m, m, 1.0_dp, fm_work, fm_K_eigvec, 0.0_dp, fm_K_sqrt)
1102 :
1103 : ! H_mp = sum_no S_mn Mr_no S_op with Mr_no <- (Mr_no + Mr_on)/2, then T holds its eigenvectors
1104 34 : CALL cp_fm_to_fm_submat(fm_red_M, fm_Mr, m, m, 1, 1, 1, 1)
1105 34 : CALL symmetrise_in_place(fm_Mr, fm_H_eigvec)
1106 34 : CALL parallel_gemm('N', 'N', m, m, m, 1.0_dp, fm_Mr, fm_K_sqrt, 0.0_dp, fm_work)
1107 34 : CALL parallel_gemm('N', 'N', m, m, m, 1.0_dp, fm_K_sqrt, fm_work, 0.0_dp, fm_H)
1108 34 : lambda(:) = 0.0_dp
1109 34 : CALL choose_eigv_solver(fm_H, fm_H_eigvec, lambda)
1110 : ! H is congruent to Mr = sum_ia b (M b), so λ_1 <= 0 means (A+B) indefinite on the span of b, or
1111 : ! round-off in a nearly dependent basis, which extension_deviation excludes at every extension
1112 34 : IF (lambda(1) <= 0.0_dp) THEN
1113 : CALL cp_abort(__LOCATION__, &
1114 : "BSE Davidson (OLSEN_STRATMANN): negative squared excitation energy. Matrix "// &
1115 : "(A+B) is not positive definite on the trial space, or the basis lost its "// &
1116 0 : "orthonormality.")
1117 : END IF
1118 272 : omega(1:n_act) = SQRT(lambda(1:n_act))
1119 :
1120 : ! R_mk = sum_n S_mn T_nk/√Ω_k, L_mk = sum_n Mr_mn R_nk/Ω_k
1121 34 : NULLIFY (fm_struct)
1122 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1123 34 : nrow_global=m, ncol_global=n_act)
1124 34 : CALL cp_fm_create(fm_R, fm_struct, name="bse_paired_R")
1125 34 : CALL cp_fm_create(fm_L, fm_struct, name="bse_paired_L")
1126 34 : CALL cp_fm_struct_release(fm_struct)
1127 34 : CALL parallel_gemm('N', 'N', m, n_act, m, 1.0_dp, fm_K_sqrt, fm_H_eigvec, 0.0_dp, fm_R)
1128 272 : CALL cp_fm_column_scale(fm_R, 1.0_dp/SQRT(omega(1:n_act)))
1129 34 : CALL parallel_gemm('N', 'N', m, n_act, m, 1.0_dp, fm_Mr, fm_R, 0.0_dp, fm_L)
1130 272 : CALL cp_fm_column_scale(fm_L, 1.0_dp/omega(1:n_act))
1131 :
1132 6936 : coef_right(:, :) = 0.0_dp
1133 6936 : coef_left(:, :) = 0.0_dp
1134 34 : CALL cp_fm_get_submatrix(fm_R, coef_right(1:m, 1:n_act))
1135 34 : CALL cp_fm_get_submatrix(fm_L, coef_left(1:m, 1:n_act))
1136 34 : CALL cp_fm_release(fm_R)
1137 34 : CALL cp_fm_release(fm_L)
1138 : END IF
1139 :
1140 34 : CALL cp_fm_release(fm_K_eigvec)
1141 34 : CALL cp_fm_release(fm_K_sqrt)
1142 34 : CALL cp_fm_release(fm_Mr)
1143 34 : CALL cp_fm_release(fm_H)
1144 34 : CALL cp_fm_release(fm_H_eigvec)
1145 34 : CALL cp_fm_release(fm_work)
1146 34 : DEALLOCATE (kappa, lambda)
1147 34 : CALL timestop(handle)
1148 :
1149 68 : END SUBROUTINE solve_reduced_paired
1150 :
1151 : ! **************************************************************************************************
1152 : !> \brief Makes the nt columns T that follow column m of b orthonormal to b and to each other,
1153 : !> T_ia,k <- T_ia,k - sum_m b_ia,m sum_jb b_jb,m T_jb,k, then T_ia,k <- sum_l T_ia,l (R^-1)_lk with
1154 : !> sum_p R_pk R_pl = sum_ia T_ia,k T_ia,l (cholqr2), and stores the products
1155 : !> (M b)_ia,k = sum_jb (A+B)_ia,jb b_jb,k, (K b)_ia,k = sum_jb (A-B)_ia,jb b_jb,k of the kept
1156 : !> columns from one application of A and B
1157 : !> \param mv_env the matrix-free A and B
1158 : !> \param fm_b basis; columns 1..m are kept, m+1..m+nt hold T on entry and the kept columns on exit
1159 : !> \param fm_Mb (A+B) b, extended in the same columns
1160 : !> \param fm_Kb (A-B) b, extended in the same columns
1161 : !> \param fm_scratch scratch of at least nt columns
1162 : !> \param m columns of the basis before the extension
1163 : !> \param nt on exit the number of columns kept
1164 : !> \param n_kernel counter of the kernel applications
1165 : !> \param n_dependent counter of the correction vectors dropped as dependent
1166 : ! **************************************************************************************************
1167 34 : SUBROUTINE extend_basis_os(mv_env, fm_b, fm_Mb, fm_Kb, fm_scratch, m, nt, n_kernel, n_dependent)
1168 :
1169 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1170 : TYPE(cp_fm_type), INTENT(IN) :: fm_b, fm_Mb, fm_Kb, fm_scratch
1171 : INTEGER, INTENT(IN) :: m
1172 : INTEGER, INTENT(INOUT) :: nt, n_kernel, n_dependent
1173 :
1174 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extend_basis_os'
1175 :
1176 : INTEGER :: handle
1177 :
1178 34 : CALL timeset(routineN, handle)
1179 :
1180 34 : CALL extend_orthonormal(fm_b, m, nt, mv_env%para_env, n_dependent)
1181 :
1182 34 : IF (nt > 0) THEN
1183 : ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
1184 34 : IF (extension_deviation(fm_b, fm_b, m, nt, mv_env%para_env) > ortho_tol) THEN
1185 0 : CPABORT("BSE Davidson: the basis lost orthonormality")
1186 : END IF
1187 : ! A b into Mb and B b into the scratch, then K b = A b - B b and M b = A b + B b
1188 34 : CALL bse_matvec_apply(mv_env, fm_b, m + 1, nt, fm_Mb, fm_scratch, 1)
1189 34 : CALL cp_fm_to_fm(fm_Mb, fm_Kb, nt, m + 1, m + 1)
1190 34 : CALL columns_axpy(-1.0_dp, fm_scratch, 1, fm_Kb, m + 1, nt)
1191 34 : CALL columns_axpy(1.0_dp, fm_scratch, 1, fm_Mb, m + 1, nt)
1192 34 : n_kernel = n_kernel + nt
1193 : END IF
1194 :
1195 34 : CALL timestop(handle)
1196 :
1197 34 : END SUBROUTINE extend_basis_os
1198 :
1199 : ! **************************************************************************************************
1200 : !> \brief One convergence step of a driver: dE_k = |E_k - E_k^prev|, the flags from CONVERGENCE_CRITERION,
1201 : !> the last required root, every root converged once the basis spans all transitions, the roots
1202 : !> in need, and the line of the iteration table
1203 : !> \param mp2_env CONVERGENCE_CRITERION, EPS_ENERGY and EPS_RES
1204 : !> \param energies tracked energies in Hartree, ascending; the first SIZE(res) are read
1205 : !> \param energies_prev the same of the previous iteration
1206 : !> \param res residual norms of the tracked roots
1207 : !> \param m subspace dimension
1208 : !> \param n_ov number of transitions
1209 : !> \param n_want requested states
1210 : !> \param iter iteration count
1211 : !> \param t_iter wall time at the start of the iteration
1212 : !> \param unit_nr output unit, positive on the writing rank only
1213 : !> \param dE energy changes in Hartree
1214 : !> \param conv convergence flags
1215 : !> \param n_req last required root
1216 : !> \param need roots in need
1217 : ! **************************************************************************************************
1218 100 : SUBROUTINE convergence_step(mp2_env, energies, energies_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
1219 100 : dE, conv, n_req, need)
1220 :
1221 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1222 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies, energies_prev, res
1223 : INTEGER, INTENT(IN) :: m, n_ov, n_want, iter
1224 : REAL(KIND=dp), INTENT(IN) :: t_iter
1225 : INTEGER, INTENT(IN) :: unit_nr
1226 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: dE
1227 : LOGICAL, DIMENSION(:), INTENT(OUT) :: conv
1228 : INTEGER, INTENT(OUT) :: n_req
1229 : LOGICAL, DIMENSION(:), INTENT(OUT) :: need
1230 :
1231 : INTEGER :: n_act
1232 :
1233 50 : n_act = SIZE(res)
1234 400 : dE(:) = ABS(energies(1:n_act) - energies_prev(1:n_act))
1235 50 : CALL convergence_flags(mp2_env, dE, res, conv)
1236 50 : n_req = multiplet_end(energies, n_want, n_act, deg_thresh)
1237 : ! a subspace that spans all transitions is exact
1238 50 : IF (m == n_ov) conv(:) = .TRUE.
1239 50 : CALL roots_in_need(conv, energies, res, n_req, need)
1240 : ! the masked MAXVAL is -HUGE when no root is in need; the unmasked term covers that case
1241 : CALL print_iteration(iter, m, COUNT(conv(1:n_want)), &
1242 : MAX(MAXVAL(res(1:n_req)), MAXVAL(res, MASK=need)), &
1243 856 : MAXVAL(dE(1:n_req)), energies(1), m_walltime() - t_iter, unit_nr)
1244 :
1245 50 : END SUBROUTINE convergence_step
1246 :
1247 : ! **************************************************************************************************
1248 : !> \brief Convergence flag of every tracked root k from CONVERGENCE_CRITERION: EN dE_k < EPS_ENERGY,
1249 : !> RES |r_k| < EPS_RES, OR either of the two, AND both
1250 : !> \param mp2_env CONVERGENCE_CRITERION, EPS_ENERGY and EPS_RES
1251 : !> \param dE energy changes in Hartree
1252 : !> \param res residual norms
1253 : !> \param conv one flag per tracked root
1254 : ! **************************************************************************************************
1255 50 : SUBROUTINE convergence_flags(mp2_env, dE, res, conv)
1256 :
1257 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1258 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dE, res
1259 : LOGICAL, DIMENSION(:), INTENT(OUT) :: conv
1260 :
1261 : INTEGER :: k
1262 :
1263 400 : DO k = 1, SIZE(conv)
1264 50 : SELECT CASE (mp2_env%bse%convergence_criterion)
1265 : CASE (bse_iter_en_cond)
1266 0 : conv(k) = dE(k) < mp2_env%bse%eps_energy
1267 : CASE (bse_iter_res_cond)
1268 0 : conv(k) = res(k) < mp2_env%bse%eps_res
1269 : CASE (bse_iter_or_cond)
1270 0 : conv(k) = dE(k) < mp2_env%bse%eps_energy .OR. res(k) < mp2_env%bse%eps_res
1271 : CASE (bse_iter_and_cond)
1272 650 : conv(k) = dE(k) < mp2_env%bse%eps_energy .AND. res(k) < mp2_env%bse%eps_res
1273 : END SELECT
1274 : END DO
1275 :
1276 50 : END SUBROUTINE convergence_flags
1277 :
1278 : ! **************************************************************************************************
1279 : !> \brief Last member of the multiplet that contains state n_first: the scan stops at the first gap of
1280 : !> at least thresh, or at n_limit. A multiplet cut by NUM_EXC_EN is thereby taken as a whole
1281 : !> \param energies ascending energies in Hartree
1282 : !> \param n_first first state of the scan
1283 : !> \param n_limit last state the scan may reach
1284 : !> \param thresh gap that ends a multiplet, in Hartree
1285 : !> \return the last member, n_first <= n_last <= n_limit
1286 : ! **************************************************************************************************
1287 56 : PURE FUNCTION multiplet_end(energies, n_first, n_limit, thresh) RESULT(n_last)
1288 :
1289 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies
1290 : INTEGER, INTENT(IN) :: n_first, n_limit
1291 : REAL(KIND=dp), INTENT(IN) :: thresh
1292 : INTEGER :: n_last
1293 :
1294 56 : n_last = n_first
1295 56 : DO WHILE (n_last < n_limit)
1296 56 : IF (energies(n_last + 1) - energies(n_last) >= thresh) EXIT
1297 56 : n_last = n_last + 1
1298 : END DO
1299 :
1300 56 : END FUNCTION multiplet_end
1301 : ! **************************************************************************************************
1302 : !> \brief Roots that still need correction vectors: the unconverged ones up to n_req, and a tracked
1303 : !> root k above them as long as E_k - |r_k| < E_n_req. An eigenvalue of a symmetric matrix
1304 : !> lies within |r_k| of E_k, so such a root can still belong to the requested states; the ABBA
1305 : !> solvers, whose residuals are not those of a symmetric problem, use the same rule as a heuristic
1306 : !> \param conv convergence flags of the tracked roots
1307 : !> \param energies tracked energies in Hartree, ascending
1308 : !> \param res residual norms
1309 : !> \param n_req last required root
1310 : !> \param need one flag per tracked root, set here
1311 : ! **************************************************************************************************
1312 50 : SUBROUTINE roots_in_need(conv, energies, res, n_req, need)
1313 :
1314 : LOGICAL, DIMENSION(:), INTENT(IN) :: conv
1315 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies, res
1316 : INTEGER, INTENT(IN) :: n_req
1317 : LOGICAL, DIMENSION(:), INTENT(OUT) :: need
1318 :
1319 : INTEGER :: k
1320 :
1321 400 : DO k = 1, SIZE(need)
1322 350 : need(k) = .NOT. conv(k)
1323 590 : IF (k > n_req) need(k) = need(k) .AND. energies(k) - res(k) < energies(n_req)
1324 : END DO
1325 :
1326 50 : END SUBROUTINE roots_in_need
1327 :
1328 : ! **************************************************************************************************
1329 : !> \brief The lowest roots in need, which receive a correction vector
1330 : !> \param mp2_env CONVERGENCE_CRITERION and EPS_RES
1331 : !> \param need roots in need, from roots_in_need
1332 : !> \param res residual norms
1333 : !> \param n_max largest number of roots
1334 : !> \param selected indices of the selected roots in ascending order, zero beyond nt
1335 : !> \param nt number of selected roots
1336 : ! **************************************************************************************************
1337 44 : SUBROUTINE select_roots(mp2_env, need, res, n_max, selected, nt)
1338 :
1339 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1340 : LOGICAL, DIMENSION(:), INTENT(IN) :: need
1341 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: res
1342 : INTEGER, INTENT(IN) :: n_max
1343 : INTEGER, DIMENSION(:), INTENT(OUT) :: selected
1344 : INTEGER, INTENT(OUT) :: nt
1345 :
1346 : INTEGER :: k
1347 :
1348 352 : selected(:) = 0
1349 44 : nt = 0
1350 352 : DO k = 1, SIZE(need)
1351 308 : IF (.NOT. need(k)) CYCLE
1352 : ! a root below EPS_RES only waits for its energy criterion
1353 110 : IF (mp2_env%bse%convergence_criterion == bse_iter_and_cond .AND. &
1354 : res(k) < mp2_env%bse%eps_res) CYCLE
1355 110 : IF (nt == n_max) EXIT
1356 110 : nt = nt + 1
1357 352 : selected(nt) = k
1358 : END DO
1359 :
1360 44 : END SUBROUTINE select_roots
1361 :
1362 : ! **************************************************************************************************
1363 : !> \brief Title and column heads of the iteration table
1364 : !> \param title line printed above the column heads
1365 : !> \param unit_nr output unit, positive on the writing rank only
1366 : ! **************************************************************************************************
1367 6 : SUBROUTINE print_iteration_header(title, unit_nr)
1368 :
1369 : CHARACTER(LEN=*), INTENT(IN) :: title
1370 : INTEGER, INTENT(IN) :: unit_nr
1371 :
1372 6 : IF (unit_nr > 0) THEN
1373 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1374 3 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', title
1375 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1376 3 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
1377 6 : '|r| in a.u., energies in eV, t the wall time of the iteration'
1378 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1379 3 : WRITE (unit_nr, '(T2,A4,T7,A5,T14,A8,T24,A6,T33,A12,T47,A12,T61,A12,T74,A7)') 'BSE|', &
1380 6 : 'Iter.', 'Z-space', 'Conv.', 'Max |r|', 'Max dE', 'Lowest E', 't (s)'
1381 : END IF
1382 :
1383 6 : END SUBROUTINE print_iteration_header
1384 :
1385 : ! **************************************************************************************************
1386 : !> \brief One line of the iteration table
1387 : !> \param iter iteration count
1388 : !> \param m subspace dimension
1389 : !> \param n_conv converged roots among the requested ones
1390 : !> \param max_res largest residual norm over the required roots and the roots in need
1391 : !> \param max_dE in Hartree
1392 : !> \param lowest_E in Hartree
1393 : !> \param t_iter wall time of the iteration in seconds
1394 : !> \param unit_nr output unit, positive on the writing rank only
1395 : ! **************************************************************************************************
1396 50 : SUBROUTINE print_iteration(iter, m, n_conv, max_res, max_dE, lowest_E, t_iter, unit_nr)
1397 :
1398 : INTEGER, INTENT(IN) :: iter, m, n_conv
1399 : REAL(KIND=dp), INTENT(IN) :: max_res, max_dE, lowest_E, t_iter
1400 : INTEGER, INTENT(IN) :: unit_nr
1401 :
1402 50 : IF (unit_nr > 0) THEN
1403 : WRITE (unit_nr, '(T2,A4,T7,I5,T14,I8,T24,I6,T33,ES12.4,T47,ES12.4,T61,F12.6,T74,F7.2)') &
1404 25 : 'BSE|', iter, m, n_conv, max_res, max_dE*evolt, lowest_E*evolt, t_iter
1405 : END IF
1406 :
1407 50 : END SUBROUTINE print_iteration
1408 :
1409 : ! **************************************************************************************************
1410 : !> \brief Table of the requested roots and of the tracked roots in need, then the abort for MAX_ITER
1411 : !> \param energies in Hartree
1412 : !> \param res residual norms
1413 : !> \param conv convergence flags
1414 : !> \param need roots in need
1415 : !> \param n_req last required root
1416 : !> \param unit_nr output unit, positive on the writing rank only
1417 : ! **************************************************************************************************
1418 0 : SUBROUTINE abort_unconverged(energies, res, conv, need, n_req, unit_nr)
1419 :
1420 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies, res
1421 : LOGICAL, DIMENSION(:), INTENT(IN) :: conv, need
1422 : INTEGER, INTENT(IN) :: n_req, unit_nr
1423 :
1424 : INTEGER :: k
1425 :
1426 0 : IF (unit_nr > 0) THEN
1427 0 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1428 0 : WRITE (unit_nr, '(T2,A4,T7,A12,T30,A11,T50,A10,T67,A14)') 'BSE|', &
1429 0 : 'Excitation n', 'Energy (eV)', '|r| (a.u.)', 'Converged'
1430 0 : DO k = 1, SIZE(need)
1431 0 : IF (k > n_req .AND. .NOT. need(k)) CYCLE
1432 0 : WRITE (unit_nr, '(T2,A4,T7,I12,T27,F14.6,T46,ES14.4,T67,L14)') 'BSE|', &
1433 0 : k, energies(k)*evolt, res(k), conv(k)
1434 : END DO
1435 : END IF
1436 : ! unconverged roots would enter postprocess_bse as converged states; the table above is the diagnostic
1437 0 : CPABORT("BSE Davidson: MAX_ITER reached before convergence")
1438 :
1439 0 : END SUBROUTINE abort_unconverged
1440 :
1441 : ! **************************************************************************************************
1442 : !> \brief Debug table of all tracked roots of a converged solve
1443 : !> \param energies in Hartree
1444 : !> \param res residual norms
1445 : !> \param n_want requested states
1446 : !> \param n_act tracked states
1447 : !> \param unit_nr output unit, positive on the writing rank only
1448 : ! **************************************************************************************************
1449 0 : SUBROUTINE print_tracked_roots(energies, res, n_want, n_act, unit_nr)
1450 :
1451 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies, res
1452 : INTEGER, INTENT(IN) :: n_want, n_act, unit_nr
1453 :
1454 : INTEGER :: k
1455 :
1456 0 : IF (unit_nr > 0) THEN
1457 0 : WRITE (unit_nr, '(T2,A10,T13,A13,T27,A9,T45,A18,T67,A14)') 'BSE|DEBUG|', &
1458 0 : 'Tracked state', 'Requested', 'Energy (eV)', '|r| (a.u.)'
1459 0 : DO k = 1, n_act
1460 0 : WRITE (unit_nr, '(T2,A10,T13,I13,T27,L9,T45,F18.10,T67,ES14.4)') 'BSE|DEBUG|', &
1461 0 : k, k <= n_want, energies(k)*evolt, res(k)
1462 : END DO
1463 : END IF
1464 :
1465 0 : END SUBROUTINE print_tracked_roots
1466 :
1467 : ! **************************************************************************************************
1468 : !> \brief Counters of a converged solve, the multiplet note and the warning for a degenerate group
1469 : !> that reaches the last buffer state
1470 : !> \param iter iterations run
1471 : !> \param n_kernel kernel applications
1472 : !> \param n_restart thick restarts
1473 : !> \param n_want requested states
1474 : !> \param n_req last required root
1475 : !> \param n_act tracked states
1476 : !> \param n_ov number of transitions
1477 : !> \param unit_nr output unit, positive on the writing rank only
1478 : !> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
1479 : !> \param n_dependent correction vectors dropped as dependent, printed when not zero
1480 : ! **************************************************************************************************
1481 6 : SUBROUTINE print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, &
1482 : n_dependent)
1483 :
1484 : INTEGER, INTENT(IN) :: iter, n_kernel, n_restart, n_want, &
1485 : n_req, n_act, n_ov, unit_nr
1486 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: ab_margin
1487 : INTEGER, INTENT(IN), OPTIONAL :: n_dependent
1488 :
1489 6 : IF (unit_nr > 0) THEN
1490 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1491 3 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Davidson iterations', iter
1492 3 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Kernel applications', n_kernel
1493 3 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Thick restarts', n_restart
1494 3 : IF (PRESENT(n_dependent)) THEN
1495 3 : IF (n_dependent > 0) THEN
1496 0 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Correction vectors dropped as dependent', n_dependent
1497 : END IF
1498 : END IF
1499 3 : IF (PRESENT(ab_margin)) THEN
1500 2 : WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
1501 4 : 'A-B margin on the trial space (eV)', ab_margin*evolt
1502 : END IF
1503 3 : IF (n_req > n_want) THEN
1504 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
1505 0 : 'Note: the requested states end inside a degenerate multiplet.'
1506 : END IF
1507 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1508 : END IF
1509 :
1510 6 : IF (n_req > n_want .AND. n_req == n_act .AND. n_act < n_ov .AND. unit_nr > 0) THEN
1511 : CALL cp_warn(__LOCATION__, &
1512 : "BSE Davidson: the degenerate group at NUM_EXC_EN reaches the last buffer state, "// &
1513 0 : "further members may lie above it. Increase BSE_ITERAT%NUM_BUFFER_STATES.")
1514 : END IF
1515 :
1516 6 : END SUBROUTINE print_summary
1517 :
1518 : ! **************************************************************************************************
1519 : !> \brief Number of wanted and tracked states and the block size, clamped to the problem size
1520 : !> \param mp2_env NUM_EXC_EN, NUM_BUFFER_STATES, BLOCK_SIZE and MAX_ITER
1521 : !> \param n_ov number of transitions
1522 : !> \param unit_nr output unit, positive on the writing rank only
1523 : !> \param n_want NUM_EXC_EN, at most n_ov
1524 : !> \param n_act wanted plus buffer states
1525 : !> \param block_size BLOCK_SIZE, or MIN(32, n_act) for -1; at most n_act
1526 : ! **************************************************************************************************
1527 6 : SUBROUTINE davidson_sizes(mp2_env, n_ov, unit_nr, n_want, n_act, block_size)
1528 :
1529 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1530 : INTEGER, INTENT(IN) :: n_ov, unit_nr
1531 : INTEGER, INTENT(OUT) :: n_want, n_act, block_size
1532 :
1533 6 : IF (mp2_env%bse%num_exc_en < 1) CPABORT("BSE_ITERAT%NUM_EXC_EN must be at least 1")
1534 6 : IF (mp2_env%bse%num_buffer_states < 0) CPABORT("BSE_ITERAT%NUM_BUFFER_STATES must not be negative")
1535 6 : IF (mp2_env%bse%block_size < 1 .AND. mp2_env%bse%block_size /= -1) THEN
1536 0 : CPABORT("BSE_ITERAT%BLOCK_SIZE must be at least 1, or -1 for the default")
1537 : END IF
1538 6 : IF (mp2_env%bse%max_iter < 1) CPABORT("BSE_ITERAT%MAX_ITER must be at least 1")
1539 :
1540 6 : n_want = mp2_env%bse%num_exc_en
1541 6 : IF (n_want > n_ov) THEN
1542 0 : IF (unit_nr > 0) CALL cp_warn(__LOCATION__, &
1543 0 : "BSE_ITERAT%NUM_EXC_EN exceeds the number of transitions and is reduced to it.")
1544 0 : n_want = n_ov
1545 : END IF
1546 6 : n_act = MIN(n_want + mp2_env%bse%num_buffer_states, n_ov)
1547 6 : block_size = mp2_env%bse%block_size
1548 : ! 32 was within 5 % of the best width on the L8 deck at 200 states; wider passes cost restarts
1549 6 : IF (block_size == -1) block_size = MIN(32, n_act)
1550 6 : IF (block_size > n_act) THEN
1551 0 : IF (unit_nr > 0) CALL cp_warn(__LOCATION__, &
1552 0 : "BSE_ITERAT%BLOCK_SIZE exceeds the number of tracked states and is reduced to it.")
1553 0 : block_size = n_act
1554 : END IF
1555 :
1556 6 : END SUBROUTINE davidson_sizes
1557 :
1558 : ! **************************************************************************************************
1559 : !> \brief Largest subspace dimension m_max: MAX_SUBSPACE_FACTOR x n_act, or 20 x n_act for -1, with
1560 : !> the solver's arrays at m_max checked against the memory budget per rank as MEMORY_CHECK says;
1561 : !> CLAMP lowers m_max to the largest that fits, never below 3 x n_act (4 x n_act for OLSEN_STRATMANN)
1562 : !> \param mp2_env MAX_SUBSPACE_FACTOR, MEMORY_CHECK and MEMORY_BUDGET_GB
1563 : !> \param mv_env n_ov, the RI share per rank and block_cols
1564 : !> \param driver driver_tda, driver_mk or driver_os
1565 : !> \param n_act tracked states
1566 : !> \param block_size correction vectors per iteration
1567 : !> \param unit_nr output unit, positive on the writing rank only
1568 : !> \param m_max the ceiling, 3 or 4 x n_act <= m_max <= n_ov
1569 : ! **************************************************************************************************
1570 12 : SUBROUTINE subspace_ceiling(mp2_env, mv_env, driver, n_act, block_size, unit_nr, m_max)
1571 :
1572 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1573 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1574 : INTEGER, INTENT(IN) :: driver, n_act, block_size, unit_nr
1575 : INTEGER, INTENT(OUT) :: m_max
1576 :
1577 : CHARACTER(LEN=16) :: bud_str, fac_str, mem_str
1578 : INTEGER :: m_floor, m_hi, m_lo, m_mid, m_start, &
1579 : min_fac, n_ri_max, nb
1580 : LOGICAL :: from_free, over, skipped
1581 : REAL(KIND=dp) :: budget_GB, dist_GB, mem_avail_GB, repl_GB
1582 :
1583 : ! the smallest ceiling that leaves room for the restart and one block: TDA and MK restart onto at
1584 : ! most 2 n_act vectors and add a block of at most n_act; OS restarts onto right and left pairs
1585 : ! under the cap m_max - 2 block_size, which has to hold the current 2 n_act
1586 6 : IF (driver == driver_os) THEN
1587 2 : min_fac = 4
1588 : ELSE
1589 4 : min_fac = 3
1590 : END IF
1591 6 : WRITE (fac_str, '(I16)') min_fac
1592 :
1593 : ! the start value: the keyword, or for -1 the generous 20 x n_act that the memory check below trims to what fits
1594 6 : IF (mp2_env%bse%max_subspace_factor == -1) THEN
1595 0 : m_start = 20*n_act
1596 6 : ELSE IF (mp2_env%bse%max_subspace_factor >= min_fac) THEN
1597 6 : m_start = mp2_env%bse%max_subspace_factor*n_act
1598 : ELSE
1599 : CALL cp_abort(__LOCATION__, &
1600 : "BSE_ITERAT%MAX_SUBSPACE_FACTOR must be at least "//TRIM(ADJUSTL(fac_str))// &
1601 0 : " (or -1) for the chosen solver")
1602 : END IF
1603 6 : m_floor = MIN(min_fac*n_act, mv_env%n_ov)
1604 6 : m_start = MIN(MAX(m_start, m_floor), mv_env%n_ov)
1605 6 : m_max = m_start
1606 :
1607 : ! the widest kernel call of the driver caps the pass width; the largest RI share of any rank
1608 : ! keeps the estimate, and so m_max, identical on all ranks
1609 2 : SELECT CASE (driver)
1610 : CASE (driver_tda)
1611 2 : nb = n_act
1612 : CASE (driver_mk)
1613 2 : nb = block_size
1614 : CASE DEFAULT
1615 6 : nb = 2*block_size
1616 : END SELECT
1617 6 : IF (mv_env%block_cols > 0) nb = MIN(nb, mv_env%block_cols)
1618 6 : n_ri_max = mv_env%n_ri_loc
1619 6 : CALL mv_env%para_env%max(n_ri_max)
1620 :
1621 : ! the budget: as given, or the share mem_fraction of the free memory per rank
1622 6 : from_free = mp2_env%bse%memory_budget_gb < 0.0_dp
1623 6 : skipped = .FALSE.
1624 6 : budget_GB = mp2_env%bse%memory_budget_gb
1625 6 : IF (from_free) THEN
1626 6 : CALL mp_mem_avail_per_rank_GB(mv_env%para_env, mem_avail_GB)
1627 6 : budget_GB = mem_fraction*mem_avail_GB
1628 6 : skipped = mem_avail_GB <= 0.0_dp
1629 6 : IF (skipped .AND. mp2_env%bse%memory_check /= bse_memcheck_off .AND. unit_nr > 0) THEN
1630 0 : CALL cp_warn(__LOCATION__, "BSE Davidson: free memory not detectable, memory check skipped")
1631 : END IF
1632 : END IF
1633 :
1634 6 : CALL davidson_footprint(driver, mv_env, m_max, n_act, block_size, nb, n_ri_max, dist_GB, repl_GB)
1635 : over = .NOT. skipped .AND. mp2_env%bse%memory_check /= bse_memcheck_off .AND. &
1636 6 : dist_GB + repl_GB > budget_GB
1637 :
1638 0 : IF (over .AND. mp2_env%bse%memory_check == bse_memcheck_clamp) THEN
1639 0 : CALL davidson_footprint(driver, mv_env, m_floor, n_act, block_size, nb, n_ri_max, dist_GB, repl_GB)
1640 0 : IF (dist_GB + repl_GB > budget_GB) THEN
1641 0 : WRITE (mem_str, '(F16.3)') dist_GB + repl_GB
1642 0 : WRITE (bud_str, '(F16.3)') budget_GB
1643 : CALL cp_abort(__LOCATION__, &
1644 : "BSE Davidson: not enough memory for the smallest subspace of "//TRIM(ADJUSTL(fac_str))// &
1645 : " vectors per state: "//TRIM(ADJUSTL(mem_str))//" GB per MPI rank, the budget is "// &
1646 0 : TRIM(ADJUSTL(bud_str))//" GB. Raise MEMORY_BUDGET_GB or use more MPI ranks.")
1647 : END IF
1648 : ! the footprint grows with m: bisect for the largest m within the budget
1649 : m_lo = m_floor
1650 : m_hi = m_start
1651 0 : DO WHILE (m_hi - m_lo > 1)
1652 0 : m_mid = (m_lo + m_hi)/2
1653 0 : CALL davidson_footprint(driver, mv_env, m_mid, n_act, block_size, nb, n_ri_max, dist_GB, repl_GB)
1654 0 : IF (dist_GB + repl_GB <= budget_GB) THEN
1655 : m_lo = m_mid
1656 : ELSE
1657 0 : m_hi = m_mid
1658 : END IF
1659 : END DO
1660 0 : m_max = m_lo
1661 0 : CALL davidson_footprint(driver, mv_env, m_max, n_act, block_size, nb, n_ri_max, dist_GB, repl_GB)
1662 : END IF
1663 :
1664 6 : IF (unit_nr > 0) THEN
1665 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1666 3 : IF (mp2_env%bse%max_subspace_factor == -1) THEN
1667 0 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension (20 x buffered states)', &
1668 0 : m_start
1669 : ELSE
1670 3 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension (MAX_SUBSPACE_FACTOR)', &
1671 6 : m_start
1672 : END IF
1673 3 : IF (m_max < m_start) THEN
1674 0 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension after the memory check', &
1675 0 : m_max
1676 : END IF
1677 3 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory per MPI rank at the subspace ceiling (GB)'
1678 3 : WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Distributed arrays', dist_GB
1679 3 : WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Replicated arrays', repl_GB
1680 3 : WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Total', dist_GB + repl_GB
1681 3 : IF (.NOT. skipped) THEN
1682 3 : IF (from_free) THEN
1683 3 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Fraction of the free memory made available', &
1684 6 : mem_fraction
1685 : END IF
1686 3 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Memory budget per MPI rank (GB)', budget_GB
1687 : END IF
1688 3 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1689 3 : IF (mp2_env%bse%memory_check == bse_memcheck_off) THEN
1690 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: OFF'
1691 3 : ELSE IF (skipped) THEN
1692 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: skipped, free memory not detectable'
1693 : ELSE
1694 0 : SELECT CASE (mp2_env%bse%memory_check)
1695 : CASE (bse_memcheck_warn)
1696 0 : IF (over) THEN
1697 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: WARN, the total exceeds the budget'
1698 : ELSE
1699 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no WARN raised'
1700 : END IF
1701 : CASE (bse_memcheck_clamp)
1702 3 : IF (over) THEN
1703 0 : WRITE (unit_nr, '(T2,A4,T7,A,I0,A,I0)') 'BSE|', &
1704 0 : 'Memory check: CLAMP applied, subspace ceiling reduced from ', m_start, ' to ', m_max
1705 : ELSE
1706 3 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no CLAMP applied'
1707 : END IF
1708 : CASE (bse_memcheck_abort)
1709 3 : IF (over) THEN
1710 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: ABORT, the total exceeds the budget'
1711 : ELSE
1712 0 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no ABORT raised'
1713 : END IF
1714 : END SELECT
1715 : END IF
1716 : END IF
1717 :
1718 6 : IF (over) THEN
1719 0 : WRITE (mem_str, '(F16.3)') dist_GB + repl_GB
1720 0 : WRITE (bud_str, '(F16.3)') budget_GB
1721 0 : SELECT CASE (mp2_env%bse%memory_check)
1722 : CASE (bse_memcheck_warn)
1723 0 : IF (unit_nr > 0) CALL cp_warn(__LOCATION__, &
1724 : "BSE Davidson: the solver's arrays at the subspace ceiling need "// &
1725 : TRIM(ADJUSTL(mem_str))//" GB per MPI rank, above the budget of "// &
1726 0 : TRIM(ADJUSTL(bud_str))//" GB.")
1727 : CASE (bse_memcheck_clamp)
1728 0 : IF (unit_nr > 0) CALL cp_warn(__LOCATION__, &
1729 : "BSE Davidson: subspace ceiling reduced to fit the memory budget; raise "// &
1730 0 : "MEMORY_BUDGET_GB or lower MAX_SUBSPACE_FACTOR to silence this.")
1731 : CASE (bse_memcheck_abort)
1732 : CALL cp_abort(__LOCATION__, &
1733 : "BSE Davidson: the solver's arrays at the subspace ceiling need "//TRIM(ADJUSTL(mem_str))// &
1734 : " GB per MPI rank, the budget is "//TRIM(ADJUSTL(bud_str))//" GB. Raise MEMORY_BUDGET_GB, "// &
1735 0 : "lower MAX_SUBSPACE_FACTOR, or set MEMORY_CHECK CLAMP.")
1736 : END SELECT
1737 : END IF
1738 :
1739 6 : END SUBROUTINE subspace_ceiling
1740 :
1741 : ! **************************************************************************************************
1742 : !> \brief Memory per MPI rank of the solver's arrays at subspace dimension m, split into the family
1743 : !> spread over the ranks and the family every rank holds whole. The counts are those of the
1744 : !> ALLOCATE and cp_fm_create statements of the drivers at a thick restart, their peak:
1745 : !> distributed 8 B (N_ov (n_basis m + n_work + n_scratch) + (n_red + n_solve) m^2)/P,
1746 : !> replicated 8 B (c_coef m n_act + m block_size + n_mv N_ov nb + n_ri nb + n_virt virt^2)
1747 : !> \param driver driver_tda, driver_mk or driver_os
1748 : !> \param mv_env n_ov, virt and the communicator
1749 : !> \param m subspace dimension
1750 : !> \param n_act tracked states
1751 : !> \param block_size correction vectors per iteration
1752 : !> \param nb kernel pass width
1753 : !> \param n_ri RI functions of a rank, the largest over the ranks
1754 : !> \param dist_GB distributed family, falls as 1/P
1755 : !> \param repl_GB replicated family, flat in P
1756 : ! **************************************************************************************************
1757 6 : SUBROUTINE davidson_footprint(driver, mv_env, m, n_act, block_size, nb, n_ri, dist_GB, repl_GB)
1758 :
1759 : INTEGER, INTENT(IN) :: driver
1760 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1761 : INTEGER, INTENT(IN) :: m, n_act, block_size, nb, n_ri
1762 : REAL(KIND=dp), INTENT(OUT) :: dist_GB, repl_GB
1763 :
1764 : INTEGER :: c_coef, n_basis, n_mv, n_red, n_scratch, &
1765 : n_solve, n_virt, n_work
1766 : REAL(KIND=dp) :: m_real, n_ov_real
1767 :
1768 : ! distributed 8 B (N_ov (n_basis m + n_work + n_scratch) + (n_red + n_solve) m^2)/P: the trial vector
1769 : ! arrays, the work and scratch columns, the reduced and the per-solve m x m matrices;
1770 : ! replicated 8 B (c_coef m n_act + m block_size + n_mv N_ov nb + n_ri nb + n_virt virt^2): the m x n_act
1771 : ! coefficient blocks with the restart work array, the pass buffers of the kernel, virt x virt for B
1772 8 : SELECT CASE (driver)
1773 : CASE (driver_tda)
1774 2 : n_basis = 2
1775 2 : n_work = MIN(2*n_act, m)
1776 2 : n_scratch = 0
1777 2 : n_red = 1
1778 2 : n_solve = 2
1779 2 : c_coef = 8
1780 2 : n_mv = 3
1781 2 : n_virt = 0
1782 : CASE (driver_mk)
1783 2 : n_basis = 3
1784 2 : n_work = MIN(2*n_act, m)
1785 2 : n_scratch = block_size
1786 2 : n_red = 1
1787 2 : n_solve = 2
1788 2 : c_coef = 10
1789 2 : n_mv = 4
1790 2 : n_virt = 1
1791 : CASE (driver_os)
1792 2 : n_basis = 3
1793 2 : n_work = MAX(2*n_act, MIN(4*n_act, m))
1794 2 : n_scratch = 2*block_size
1795 2 : n_red = 2
1796 2 : n_solve = 6
1797 2 : c_coef = 15
1798 2 : n_mv = 4
1799 2 : n_virt = 1
1800 : CASE DEFAULT
1801 : ! unreachable for the three driver constants; the DEFAULT keeps every counter defined
1802 6 : CPABORT("BSE Davidson: unknown driver in the memory estimate")
1803 : END SELECT
1804 :
1805 : ! in REAL: m^2 overflows INTEGER(4) above 46340
1806 6 : m_real = REAL(m, dp)
1807 6 : n_ov_real = REAL(mv_env%n_ov, dp)
1808 : dist_GB = 8.0E-9_dp*(n_ov_real*REAL(n_basis*m + n_work + n_scratch, dp) + REAL(n_red + n_solve, dp)*m_real*m_real)/ &
1809 6 : REAL(mv_env%para_env%num_pe, dp)
1810 : repl_GB = 8.0E-9_dp*(REAL(c_coef, dp)*m_real*REAL(n_act, dp) + m_real*REAL(block_size, dp) + &
1811 : REAL(n_mv, dp)*n_ov_real*REAL(nb, dp) + REAL(n_ri, dp)*REAL(nb, dp) + &
1812 6 : REAL(n_virt, dp)*REAL(mv_env%virt, dp)**2)
1813 :
1814 6 : END SUBROUTINE davidson_footprint
1815 :
1816 : ! **************************************************************************************************
1817 : !> \brief Unit vectors Z_ia,k = δ_ia,diag_order(k) on the entries diag_order(1..n_used) with the lowest diagonal, never
1818 : !> cutting a degenerate group. The sorted diagonal is replicated; each rank writes the entries
1819 : !> that fall on its own rows
1820 : !> \param diag the diagonal chosen by PRECONDITIONER, replicated
1821 : !> \param n_guess guess vectors wanted
1822 : !> \param n_max largest number of guess vectors
1823 : !> \param deg_thresh_diag entries of the sorted diagonal closer than this form one degenerate group
1824 : !> \param fm_Z receives the unit vectors in its columns 1..n_used
1825 : !> \param n_used guess vectors written, n_guess <= n_used <= n_max
1826 : !> \param n_given leading columns that the caller has filled; they replace the lowest unit vectors
1827 : ! **************************************************************************************************
1828 6 : SUBROUTINE initial_guess(diag, n_guess, n_max, deg_thresh_diag, fm_Z, n_used, n_given)
1829 :
1830 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag
1831 : INTEGER, INTENT(IN) :: n_guess, n_max
1832 : REAL(KIND=dp), INTENT(IN) :: deg_thresh_diag
1833 : TYPE(cp_fm_type), INTENT(IN) :: fm_Z
1834 : INTEGER, INTENT(OUT) :: n_used
1835 : INTEGER, INTENT(IN), OPTIONAL :: n_given
1836 :
1837 : INTEGER :: iloc, k, k_first, n_ov, nrow_local
1838 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: diag_order
1839 6 : INTEGER, DIMENSION(:), POINTER :: row_indices
1840 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_sorted
1841 :
1842 6 : n_ov = SIZE(diag)
1843 30 : ALLOCATE (diag_sorted(n_ov), diag_order(n_ov))
1844 294 : diag_sorted(:) = diag(:)
1845 6 : CALL sort(diag_sorted, n_ov, diag_order)
1846 :
1847 6 : n_used = multiplet_end(diag_sorted, n_guess, n_max, deg_thresh_diag)
1848 :
1849 6 : k_first = 1
1850 6 : IF (PRESENT(n_given)) k_first = n_given + 1
1851 :
1852 6 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
1853 48 : DO k = k_first, n_used
1854 1056 : DO iloc = 1, nrow_local
1855 1050 : IF (row_indices(iloc) == diag_order(k)) fm_Z%local_data(iloc, k) = 1.0_dp
1856 : END DO
1857 : END DO
1858 :
1859 6 : DEALLOCATE (diag_sorted, diag_order)
1860 :
1861 12 : END SUBROUTINE initial_guess
1862 :
1863 : ! **************************************************************************************************
1864 : !> \brief Initial guess from the exact A (and B) on the n_sub transitions with the lowest diagonal:
1865 : !> TDA: sum_q A_sub,pq V_qk = θ_k V_pk; ABBA: x_pk = sum_q (K_sub^-1/2)_pq t_qk with
1866 : !> sum_qrs (K_sub^1/2)_pq M_sub,qr (K_sub^1/2)_rs t_sk = θ_k^2 t_pk, K = A-B, M = A+B. The Ritz
1867 : !> values θ_k are upper bounds of the true energies (Rayleigh-Ritz for TDA, Cauchy-type
1868 : !> interlacing of Bai and Li for the linear response problem), kept for check_guess_bound. A
1869 : !> degenerate group at the upper boundary is taken as a whole. The block and its eigenvectors
1870 : !> are replicated; each rank writes the entries that fall on its own rows
1871 : !> \param mv_env the matrix-free A and B, for the exact block
1872 : !> \param mp2_env NUM_GUESS_TRANSITIONS and BSE_DEBUG_PRINT
1873 : !> \param diag the diagonal chosen by PRECONDITIONER, replicated; its lowest entries select the block
1874 : !> \param n_guess guess vectors wanted
1875 : !> \param n_max largest number of guess vectors
1876 : !> \param deg_thresh_diag Ritz values closer than this form one degenerate group
1877 : !> \param fm_Z receives the guess vectors in its columns 1..n_used
1878 : !> \param n_used guess vectors written, n_guess <= n_used <= n_max; 0 for abba_indefinite
1879 : !> \param theta_sub allocated here, the n_sub Ritz values in Hartree
1880 : !> \param unit_nr output unit, positive on the writing rank only
1881 : !> \param n_given leading columns that the caller has filled; they replace the lowest guess vectors
1882 : !> \param abba_status ABBA problem when present; abba_indefinite if K_sub is not positive definite
1883 : !> \param ab_margin smallest eigenvalue of K_sub in Hartree
1884 : !> \param n_pair ABBA: partners y_k = K_sub x_k/θ_k of the guess vectors x_k not given by the caller,
1885 : !> written to the n_pair columns after n_used; n_used + n_pair stays within n_max
1886 : ! **************************************************************************************************
1887 0 : SUBROUTINE initial_guess_subblock(mv_env, mp2_env, diag, n_guess, n_max, deg_thresh_diag, fm_Z, n_used, &
1888 : theta_sub, unit_nr, n_given, abba_status, ab_margin, n_pair)
1889 :
1890 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1891 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1892 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag
1893 : INTEGER, INTENT(IN) :: n_guess, n_max
1894 : REAL(KIND=dp), INTENT(IN) :: deg_thresh_diag
1895 : TYPE(cp_fm_type), INTENT(IN) :: fm_Z
1896 : INTEGER, INTENT(OUT) :: n_used
1897 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1898 : INTENT(OUT) :: theta_sub
1899 : INTEGER, INTENT(IN) :: unit_nr
1900 : INTEGER, INTENT(IN), OPTIONAL :: n_given
1901 : INTEGER, INTENT(OUT), OPTIONAL :: abba_status
1902 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ab_margin
1903 : INTEGER, INTENT(OUT), OPTIONAL :: n_pair
1904 :
1905 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initial_guess_subblock'
1906 :
1907 : INTEGER :: handle, iloc, k, k_first, n_ov, n_sub, &
1908 : n_top, nrow_local
1909 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: diag_order, rank_in_block
1910 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
1911 : LOGICAL :: ok
1912 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_sorted, eig_K
1913 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A_sub, B_sub, eigvec_H, eigvec_K, &
1914 0 : guess_x, guess_y, H_sub, K_mhalf, &
1915 0 : K_phalf, K_sub, M_sub
1916 : TYPE(mp_para_env_type), POINTER :: para_env
1917 :
1918 0 : CALL timeset(routineN, handle)
1919 :
1920 0 : IF (PRESENT(n_pair)) n_pair = 0
1921 0 : para_env => mv_env%para_env
1922 0 : n_ov = SIZE(diag)
1923 0 : n_sub = MIN(n_ov, MAX(mp2_env%bse%num_guess_transitions, n_guess))
1924 0 : ALLOCATE (diag_sorted(n_ov), diag_order(n_ov), rank_in_block(n_ov))
1925 0 : diag_sorted(:) = diag(:)
1926 0 : CALL sort(diag_sorted, n_ov, diag_order)
1927 :
1928 0 : ALLOCATE (A_sub(n_sub, n_sub), guess_x(n_sub, n_sub), theta_sub(n_sub))
1929 0 : ok = .TRUE.
1930 0 : IF (PRESENT(abba_status)) THEN
1931 0 : abba_status = abba_ok
1932 : ALLOCATE (B_sub(n_sub, n_sub), K_sub(n_sub, n_sub), M_sub(n_sub, n_sub), H_sub(n_sub, n_sub), &
1933 : eigvec_K(n_sub, n_sub), eigvec_H(n_sub, n_sub), K_phalf(n_sub, n_sub), &
1934 0 : K_mhalf(n_sub, n_sub), eig_K(n_sub))
1935 0 : CALL bse_matvec_subblock(mv_env, diag_order(1:n_sub), A_sub, B_sub)
1936 : ! (K_sub^±1/2)_pq = sum_r U_pr κ_r^±1/2 U_qr from K_sub,pq = sum_r U_pr κ_r U_qr
1937 0 : K_sub(:, :) = A_sub(:, :) - B_sub(:, :)
1938 0 : CALL solve_replicated(K_sub, n_sub, para_env, eig_K, eigvec_K)
1939 0 : ab_margin = eig_K(1)
1940 0 : IF (eig_K(1) <= 0.0_dp) THEN
1941 0 : abba_status = abba_indefinite
1942 0 : ok = .FALSE.
1943 0 : n_used = 0
1944 : ELSE
1945 0 : DO k = 1, n_sub
1946 0 : K_phalf(:, k) = eigvec_K(:, k)*SQRT(eig_K(k))
1947 0 : K_mhalf(:, k) = eigvec_K(:, k)/SQRT(eig_K(k))
1948 : END DO
1949 0 : K_phalf(:, :) = MATMUL(K_phalf, TRANSPOSE(eigvec_K))
1950 0 : K_mhalf(:, :) = MATMUL(K_mhalf, TRANSPOSE(eigvec_K))
1951 : ! M_sub = 2 A_sub - K_sub, H_ps = sum_qr (K^1/2)_pq M_qr (K^1/2)_rs, x_pk = sum_q (K^-1/2)_pq t_qk
1952 0 : M_sub(:, :) = 2.0_dp*A_sub(:, :) - K_sub(:, :)
1953 0 : H_sub(:, :) = MATMUL(K_phalf, MATMUL(M_sub, K_phalf))
1954 0 : CALL solve_replicated(H_sub, n_sub, para_env, theta_sub, eigvec_H)
1955 : ! the same failure as in solve_reduced_paired: with K_sub positive definite, M_sub is indefinite
1956 0 : IF (theta_sub(1) <= 0.0_dp) THEN
1957 0 : CPABORT("BSE Davidson: the guess block gives a non-positive squared excitation energy")
1958 : END IF
1959 0 : theta_sub(:) = SQRT(theta_sub(:))
1960 0 : guess_x(:, :) = MATMUL(K_mhalf, eigvec_H)
1961 0 : IF (PRESENT(n_pair)) THEN
1962 : ! partners guess_y_k = K_sub guess_x_k/θ_k
1963 0 : ALLOCATE (guess_y(n_sub, n_sub))
1964 0 : guess_y(:, :) = MATMUL(K_sub, guess_x)
1965 0 : DO k = 1, n_sub
1966 0 : guess_y(:, k) = guess_y(:, k)/theta_sub(k)
1967 : END DO
1968 : END IF
1969 : END IF
1970 0 : DEALLOCATE (B_sub, K_sub, M_sub, H_sub, eigvec_K, eigvec_H, K_phalf, K_mhalf, eig_K)
1971 : ELSE
1972 0 : CALL bse_matvec_subblock(mv_env, diag_order(1:n_sub), A_sub)
1973 0 : CALL solve_replicated(A_sub, n_sub, para_env, theta_sub, guess_x)
1974 : END IF
1975 :
1976 0 : IF (ok) THEN
1977 0 : k_first = 1
1978 0 : IF (PRESENT(n_given)) k_first = n_given + 1
1979 :
1980 0 : n_top = MIN(n_max, n_sub)
1981 : ! the partners of the columns k_first..n_used share the n_max columns with the guess
1982 0 : IF (PRESENT(n_pair)) n_top = MIN(n_top, (n_max + k_first - 1)/2)
1983 0 : n_used = multiplet_end(theta_sub, n_guess, n_top, deg_thresh_diag)
1984 :
1985 0 : rank_in_block(:) = 0
1986 0 : DO k = 1, n_sub
1987 0 : rank_in_block(diag_order(k)) = k
1988 : END DO
1989 0 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
1990 0 : DO k = k_first, n_used
1991 0 : DO iloc = 1, nrow_local
1992 0 : IF (rank_in_block(row_indices(iloc)) > 0) fm_Z%local_data(iloc, k) = guess_x(rank_in_block(row_indices(iloc)), k)
1993 : END DO
1994 : END DO
1995 0 : IF (PRESENT(n_pair)) THEN
1996 0 : n_pair = MIN(n_used - k_first + 1, n_max - n_used)
1997 0 : DO k = 1, n_pair
1998 0 : DO iloc = 1, nrow_local
1999 0 : IF (rank_in_block(row_indices(iloc)) > 0) THEN
2000 0 : fm_Z%local_data(iloc, n_used + k) = guess_y(rank_in_block(row_indices(iloc)), k_first + k - 1)
2001 : END IF
2002 : END DO
2003 : END DO
2004 : END IF
2005 :
2006 0 : IF (unit_nr > 0) THEN
2007 0 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Transitions in the exact guess block', n_sub
2008 0 : IF (mp2_env%bse%bse_debug_print) THEN
2009 0 : WRITE (unit_nr, '(T2,A10,T13,A,T67,F14.6)') 'BSE|DEBUG|', &
2010 0 : 'Lowest Ritz value of the guess block (eV)', theta_sub(1)*evolt
2011 : END IF
2012 : END IF
2013 : END IF
2014 :
2015 0 : DEALLOCATE (diag_sorted, diag_order, rank_in_block, A_sub, guess_x)
2016 0 : IF (ALLOCATED(guess_y)) DEALLOCATE (guess_y)
2017 :
2018 0 : CALL timestop(handle)
2019 :
2020 0 : END SUBROUTINE initial_guess_subblock
2021 :
2022 : ! **************************************************************************************************
2023 : !> \brief A converged energy above the Ritz value of the guess block with the same index is a
2024 : !> provable miss: both are upper bounds of the true eigenvalue with that index
2025 : !> \param theta_sub Ritz values of the guess block in Hartree, ascending
2026 : !> \param energies converged energies in Hartree
2027 : !> \param n_req roots compared, the required ones
2028 : !> \param unit_nr output unit, positive on the writing rank only
2029 : ! **************************************************************************************************
2030 0 : SUBROUTINE check_guess_bound(theta_sub, energies, n_req, unit_nr)
2031 :
2032 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: theta_sub, energies
2033 : INTEGER, INTENT(IN) :: n_req, unit_nr
2034 :
2035 : REAL(KIND=dp) :: excess
2036 :
2037 0 : excess = MAXVAL(energies(1:n_req) - theta_sub(1:n_req))
2038 0 : IF (unit_nr > 0) THEN
2039 0 : WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
2040 0 : 'Largest excess over the guess-block bound (eV)', excess*evolt
2041 : END IF
2042 0 : IF (excess > deg_thresh .AND. unit_nr > 0) THEN
2043 : CALL cp_warn(__LOCATION__, &
2044 : "A converged BSE state lies above the upper bound given by the guess block, "// &
2045 : "so a lower state was missed. Raise BSE_ITERAT%NUM_GUESS_TRANSITIONS or "// &
2046 0 : "BSE_ITERAT%NUM_BUFFER_STATES.")
2047 : END IF
2048 :
2049 0 : END SUBROUTINE check_guess_bound
2050 :
2051 : ! **************************************************************************************************
2052 : !> \brief sum_q (red_pq + red_qp)/2 c_qk = θ_k c_pk for p, q <= m, red_pq the live block of fm_red, solved
2053 : !> on the grid of blacs_env: θ_k in theta ascending, c_pk of the first n_vec k gathered into
2054 : !> coef_ritz
2055 : !> \param fm_red the reduced matrix; its live block is the first m rows and columns
2056 : !> \param m dimension of the live block
2057 : !> \param n_vec eigenvectors gathered; the drivers never read more than n_act
2058 : !> \param para_env communicator of the grid
2059 : !> \param blacs_env grid the reduced problem is solved on
2060 : !> \param theta θ_k in the first m entries, zero beyond
2061 : !> \param coef_ritz c_pk in the first m rows and n_vec columns, zero beyond
2062 : ! **************************************************************************************************
2063 40 : SUBROUTINE solve_reduced(fm_red, m, n_vec, para_env, blacs_env, theta, coef_ritz)
2064 :
2065 : TYPE(cp_fm_type), INTENT(IN) :: fm_red
2066 : INTEGER, INTENT(IN) :: m, n_vec
2067 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2068 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2069 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: theta
2070 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: coef_ritz
2071 :
2072 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_reduced'
2073 :
2074 : INTEGER :: handle, n_get
2075 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
2076 : TYPE(cp_fm_type) :: fm_reduced_eigvec, fm_reduced_sym
2077 :
2078 40 : CALL timeset(routineN, handle)
2079 :
2080 40 : CALL reduced_live_block(fm_red, m, para_env, blacs_env, "bse_reduced", fm_reduced_sym)
2081 40 : CALL cp_fm_create(fm_reduced_eigvec, fm_reduced_sym%matrix_struct, name="bse_reduced_vectors")
2082 : ! a_pq <- (red_pq + red_qp)/2
2083 40 : CALL symmetrise_in_place(fm_reduced_sym, fm_reduced_eigvec)
2084 120 : ALLOCATE (eigval(m))
2085 40 : eigval(:) = 0.0_dp
2086 40 : CALL choose_eigv_solver(fm_reduced_sym, fm_reduced_eigvec, eigval)
2087 :
2088 1006 : theta(:) = 0.0_dp
2089 606 : theta(1:m) = eigval(:)
2090 7082 : coef_ritz(:, :) = 0.0_dp
2091 40 : n_get = MIN(m, n_vec)
2092 40 : CALL cp_fm_get_submatrix(fm_reduced_eigvec, coef_ritz(1:m, 1:n_get), 1, 1, m, n_get)
2093 :
2094 40 : DEALLOCATE (eigval)
2095 40 : CALL cp_fm_release(fm_reduced_sym)
2096 40 : CALL cp_fm_release(fm_reduced_eigvec)
2097 :
2098 40 : CALL timestop(handle)
2099 :
2100 40 : END SUBROUTINE solve_reduced
2101 :
2102 : ! **************************************************************************************************
2103 : !> \brief sum_q (mat_pq + mat_qp)/2 coef_qk = θ_k coef_pk for p, q <= m, a small replicated block solved
2104 : !> on one rank and broadcast: the fixed-size blocks of the guess and of cholqr2, never the
2105 : !> reduced matrix
2106 : !> \param mat the block; its first m rows and columns are read
2107 : !> \param m dimension of the block
2108 : !> \param para_env communicator; the solve runs on its source rank
2109 : !> \param theta ascending
2110 : !> \param coef all m eigenvectors in the columns
2111 : ! **************************************************************************************************
2112 0 : SUBROUTINE solve_replicated(mat, m, para_env, theta, coef)
2113 :
2114 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: mat
2115 : INTEGER, INTENT(IN) :: m
2116 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2117 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: theta
2118 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: coef
2119 :
2120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_replicated'
2121 :
2122 : INTEGER :: handle
2123 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval
2124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_sym
2125 :
2126 0 : CALL timeset(routineN, handle)
2127 :
2128 0 : ALLOCATE (mat_sym(m, m), eigval(m))
2129 0 : mat_sym(:, :) = 0.5_dp*(mat(1:m, 1:m) + TRANSPOSE(mat(1:m, 1:m)))
2130 0 : eigval(:) = 0.0_dp
2131 0 : IF (para_env%is_source()) CALL diamat_all(mat_sym, eigval)
2132 0 : CALL para_env%bcast(mat_sym)
2133 0 : CALL para_env%bcast(eigval)
2134 :
2135 0 : theta(:) = 0.0_dp
2136 0 : theta(1:m) = eigval(:)
2137 0 : coef(:, :) = 0.0_dp
2138 0 : coef(1:m, 1:m) = mat_sym(:, :)
2139 0 : DEALLOCATE (mat_sym, eigval)
2140 :
2141 0 : CALL timestop(handle)
2142 :
2143 0 : END SUBROUTINE solve_replicated
2144 :
2145 : ! **************************************************************************************************
2146 : !> \brief a_pq <- (a_pq + a_qp)/2 through a scratch matrix of the same shape: the eigensolvers read one triangle
2147 : !> \param fm_a symmetrised in place
2148 : !> \param fm_scratch overwritten with the transpose
2149 : ! **************************************************************************************************
2150 108 : SUBROUTINE symmetrise_in_place(fm_a, fm_scratch)
2151 :
2152 : TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_scratch
2153 :
2154 108 : CALL cp_fm_transpose(fm_a, fm_scratch)
2155 108 : CALL cp_fm_scale_and_add(0.5_dp, fm_a, 0.5_dp, fm_scratch)
2156 :
2157 108 : END SUBROUTINE symmetrise_in_place
2158 :
2159 : ! **************************************************************************************************
2160 : !> \brief Orthonormal coefficient vectors of the restart basis from the candidate columns in their
2161 : !> order of priority, at most n_cap of them: per candidate k, u_p = cand_pk -
2162 : !> sum_j (sum_q coef_restart_qj cand_qk) coef_restart_pj over the vectors j kept so far, twice,
2163 : !> then coef_restart_p,new = u_p/|u| unless |u| is below coef_norm_drop
2164 : !> \param cand candidate coefficient columns, zero-padded to m rows: the current Ritz vectors, then
2165 : !> the previous ones
2166 : !> \param m rows of the coefficient vectors, the current basis size
2167 : !> \param n_cand candidate columns
2168 : !> \param n_cap largest number of vectors kept
2169 : !> \param coef_restart m x n_new restart basis
2170 : !> \param n_new vectors kept, at most MIN(n_cap, n_cand)
2171 : ! **************************************************************************************************
2172 24 : SUBROUTINE restart_basis(cand, m, n_cand, n_cap, coef_restart, n_new)
2173 :
2174 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cand
2175 : INTEGER, INTENT(IN) :: m, n_cand, n_cap
2176 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: coef_restart
2177 : INTEGER, INTENT(OUT) :: n_new
2178 :
2179 : CHARACTER(LEN=*), PARAMETER :: routineN = 'restart_basis'
2180 : REAL(KIND=dp), PARAMETER :: coef_norm_drop = 1.0E-8_dp
2181 :
2182 : INTEGER :: handle, ipass, j, k
2183 : REAL(KIND=dp) :: norm
2184 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: basis
2185 :
2186 24 : CALL timeset(routineN, handle)
2187 :
2188 : ! never more columns than candidates: the paired solver's cap is close to m_max
2189 96 : ALLOCATE (basis(m, MIN(n_cap, n_cand)))
2190 24 : n_new = 0
2191 360 : DO k = 1, n_cand
2192 352 : IF (n_new == n_cap) EXIT
2193 6440 : basis(:, n_new + 1) = cand(1:m, k)
2194 : ! a previous vector is nearly parallel to its successor, one pass leaves the basis
2195 : ! non-orthonormal at the 1e-6 level and stalls the residuals there
2196 992 : DO ipass = 1, 2
2197 4898 : DO j = 1, n_new
2198 : basis(:, n_new + 1) = basis(:, n_new + 1) - &
2199 159510 : DOT_PRODUCT(basis(:, j), basis(:, n_new + 1))*basis(:, j)
2200 : END DO
2201 12750 : norm = NORM2(basis(:, n_new + 1))
2202 664 : IF (norm < coef_norm_drop) EXIT
2203 12956 : basis(:, n_new + 1) = basis(:, n_new + 1)/norm
2204 : END DO
2205 : ! relative: a candidate that keeps less than this of its unit norm repeats the vectors kept
2206 336 : IF (norm < coef_norm_drop) CYCLE
2207 360 : n_new = n_new + 1
2208 : END DO
2209 15676 : coef_restart(:, :) = 0.0_dp
2210 6334 : coef_restart(1:m, 1:n_new) = basis(:, 1:n_new)
2211 24 : DEALLOCATE (basis)
2212 :
2213 24 : CALL timestop(handle)
2214 :
2215 24 : END SUBROUTINE restart_basis
2216 :
2217 : ! **************************************************************************************************
2218 : !> \brief out_ia,k = alpha sum_m V_ia,m coef_mk + beta out_ia,k for m <= nv, k <= nc on the local rows ia
2219 : !> On the npe x 1 grid whole columns are local, so the rotation is one local DGEMM
2220 : !> \param fm_V basis, columns 1..nv read
2221 : !> \param nv columns of fm_V read
2222 : !> \param coef nv x nc coefficients, replicated
2223 : !> \param nc columns written
2224 : !> \param fm_out receives the rotated columns
2225 : !> \param alpha scale of the product
2226 : !> \param beta scale of the previous content of fm_out
2227 : !> \param out_col first column of fm_out written, 1 by default
2228 : ! **************************************************************************************************
2229 432 : SUBROUTINE subspace_rotate(fm_V, nv, coef, nc, fm_out, alpha, beta, out_col)
2230 :
2231 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2232 : INTEGER, INTENT(IN) :: nv
2233 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2234 : INTENT(IN) :: coef
2235 : INTEGER, INTENT(IN) :: nc
2236 : TYPE(cp_fm_type), INTENT(IN) :: fm_out
2237 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
2238 : INTEGER, INTENT(IN), OPTIONAL :: out_col
2239 :
2240 : CHARACTER(LEN=*), PARAMETER :: routineN = 'subspace_rotate'
2241 :
2242 : INTEGER :: handle, nrow_local, o_col
2243 :
2244 216 : CALL timeset(routineN, handle)
2245 :
2246 216 : o_col = 1
2247 216 : IF (PRESENT(out_col)) o_col = out_col
2248 216 : CALL cp_fm_get_info(fm_V, nrow_local=nrow_local)
2249 216 : IF (nrow_local > 0) THEN
2250 : CALL DGEMM('N', 'N', nrow_local, nc, nv, alpha, fm_V%local_data, SIZE(fm_V%local_data, 1), &
2251 216 : coef, SIZE(coef, 1), beta, fm_out%local_data(:, o_col:o_col + nc - 1), SIZE(fm_out%local_data, 1))
2252 : END IF
2253 :
2254 216 : CALL timestop(handle)
2255 :
2256 216 : END SUBROUTINE subspace_rotate
2257 :
2258 : ! **************************************************************************************************
2259 : !> \brief V_ia,k <- sum_p V_ia,p coef_pk for p <= m, k <= nq through the work matrix, and the same rotation
2260 : !> of up to two further matrices that share the basis, so that a thick restart is one call
2261 : !> \param fm_V basis; on exit columns 1..nq hold the rotation and nq+1..m are zero
2262 : !> \param m columns of fm_V read
2263 : !> \param coef restart basis, m x nq
2264 : !> \param nq columns of the rotated basis
2265 : !> \param fm_work at least nq columns
2266 : !> \param fm_V2 rotated alike
2267 : !> \param fm_V3 rotated alike
2268 : ! **************************************************************************************************
2269 24 : SUBROUTINE rotate_in_place(fm_V, m, coef, nq, fm_work, fm_V2, fm_V3)
2270 :
2271 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2272 : INTEGER, INTENT(IN) :: m
2273 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coef
2274 : INTEGER, INTENT(IN) :: nq
2275 : TYPE(cp_fm_type), INTENT(IN) :: fm_work
2276 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_V2, fm_V3
2277 :
2278 24 : CALL rotate_one(fm_V)
2279 24 : IF (PRESENT(fm_V2)) CALL rotate_one(fm_V2)
2280 24 : IF (PRESENT(fm_V3)) CALL rotate_one(fm_V3)
2281 :
2282 : CONTAINS
2283 :
2284 : ! **************************************************************************************************
2285 : !> \brief The rotation of one matrix
2286 : !> \param fm rotated in place
2287 : ! **************************************************************************************************
2288 132 : SUBROUTINE rotate_one(fm)
2289 :
2290 : TYPE(cp_fm_type), INTENT(IN) :: fm
2291 :
2292 : INTEGER :: nrow_local
2293 :
2294 17726 : CALL subspace_rotate(fm, m, coef(1:m, 1:nq), nq, fm_work, 1.0_dp, 0.0_dp)
2295 66 : CALL cp_fm_to_fm(fm_work, fm, nq)
2296 66 : CALL cp_fm_get_info(fm, nrow_local=nrow_local)
2297 7716 : fm%local_data(1:nrow_local, nq + 1:m) = 0.0_dp
2298 :
2299 66 : END SUBROUTINE rotate_one
2300 :
2301 : END SUBROUTINE rotate_in_place
2302 : ! **************************************************************************************************
2303 : !> \brief T_ia,t0+k-1 <- T_ia,t0+k-1 + alpha U_ia,u0+k-1 for k <= n
2304 : !> On the npe x 1 grid whole columns are local, so the update needs no communication
2305 : !> \param alpha scale of the added columns
2306 : !> \param fm_U columns added
2307 : !> \param u0 first column read
2308 : !> \param fm_T columns updated
2309 : !> \param t0 first column written
2310 : !> \param n number of columns
2311 : ! **************************************************************************************************
2312 100 : SUBROUTINE columns_axpy(alpha, fm_U, u0, fm_T, t0, n)
2313 :
2314 : REAL(KIND=dp), INTENT(IN) :: alpha
2315 : TYPE(cp_fm_type), INTENT(IN) :: fm_U
2316 : INTEGER, INTENT(IN) :: u0
2317 : TYPE(cp_fm_type), INTENT(IN) :: fm_T
2318 : INTEGER, INTENT(IN) :: t0, n
2319 :
2320 : INTEGER :: nrow_local
2321 :
2322 100 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local)
2323 : fm_T%local_data(1:nrow_local, t0:t0 + n - 1) = fm_T%local_data(1:nrow_local, t0:t0 + n - 1) + &
2324 7800 : alpha*fm_U%local_data(1:nrow_local, u0:u0 + n - 1)
2325 :
2326 100 : END SUBROUTINE columns_axpy
2327 :
2328 : ! **************************************************************************************************
2329 : !> \brief G_kl = sum_ia U_ia,u0+k-1 V_ia,v0+l-1 for k <= nu, l <= nv, local product and sum over the ranks
2330 : !> The npe x 1 grid splits the rows over the ranks, which is what the one allreduce completes
2331 : !> \param fm_U left factor
2332 : !> \param u0 first column of fm_U
2333 : !> \param nu columns of fm_U
2334 : !> \param fm_V right factor
2335 : !> \param v0 first column of fm_V
2336 : !> \param nv columns of fm_V
2337 : !> \param para_env communicator of the row-distributed matrices
2338 : !> \param gram nu x nv, replicated
2339 : ! **************************************************************************************************
2340 840 : SUBROUTINE subspace_gram(fm_U, u0, nu, fm_V, v0, nv, para_env, gram)
2341 :
2342 : TYPE(cp_fm_type), INTENT(IN) :: fm_U
2343 : INTEGER, INTENT(IN) :: u0, nu
2344 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2345 : INTEGER, INTENT(IN) :: v0, nv
2346 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2347 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2348 : INTENT(OUT) :: gram
2349 :
2350 : CHARACTER(LEN=*), PARAMETER :: routineN = 'subspace_gram'
2351 :
2352 : INTEGER :: handle, nrow_local
2353 :
2354 420 : CALL timeset(routineN, handle)
2355 :
2356 420 : CALL cp_fm_get_info(fm_U, nrow_local=nrow_local)
2357 15744 : gram(:, :) = 0.0_dp
2358 420 : IF (nrow_local > 0) THEN
2359 : CALL DGEMM('T', 'N', nu, nv, nrow_local, 1.0_dp, &
2360 : fm_U%local_data(:, u0:u0 + nu - 1), SIZE(fm_U%local_data, 1), &
2361 : fm_V%local_data(:, v0:v0 + nv - 1), SIZE(fm_V%local_data, 1), &
2362 420 : 0.0_dp, gram, SIZE(gram, 1))
2363 : END IF
2364 420 : CALL para_env%sum(gram)
2365 :
2366 420 : CALL timestop(handle)
2367 :
2368 420 : END SUBROUTINE subspace_gram
2369 :
2370 : ! **************************************************************************************************
2371 : !> \brief Overlaps of the whole basis U with a block of nv columns of V, written into the distributed
2372 : !> reduced matrix: fm_red_k,v0+l-1 <- sum_ia U_ia,k V_ia,v0+l-1 for k <= nu, l <= nv, in column
2373 : !> blocks of nb so that the replicated temporary is nu x nb and never nu x nv
2374 : !> \param fm_U the basis, columns 1..nu
2375 : !> \param nu columns of fm_U, rows of fm_red written
2376 : !> \param fm_V the block of columns
2377 : !> \param v0 first column of fm_V read and first column of fm_red written
2378 : !> \param nv columns of fm_V read
2379 : !> \param nb width of one column block
2380 : !> \param para_env communicator of the row-distributed matrices
2381 : !> \param fm_red the distributed reduced matrix, on the diagonalisation grid
2382 : ! **************************************************************************************************
2383 10 : SUBROUTINE reduced_gram_blocks(fm_U, nu, fm_V, v0, nv, nb, para_env, fm_red)
2384 :
2385 : TYPE(cp_fm_type), INTENT(IN) :: fm_U
2386 : INTEGER, INTENT(IN) :: nu
2387 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2388 : INTEGER, INTENT(IN) :: v0, nv, nb
2389 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2390 : TYPE(cp_fm_type), INTENT(IN) :: fm_red
2391 :
2392 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reduced_gram_blocks'
2393 :
2394 : INTEGER :: c0, handle, nc
2395 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gram_block
2396 :
2397 10 : CALL timeset(routineN, handle)
2398 :
2399 40 : ALLOCATE (gram_block(nu, nb))
2400 24 : DO c0 = v0, v0 + nv - 1, nb
2401 14 : nc = MIN(nb, v0 + nv - c0)
2402 14 : CALL subspace_gram(fm_U, 1, nu, fm_V, c0, nc, para_env, gram_block(:, 1:nc))
2403 24 : CALL cp_fm_set_submatrix(fm_red, gram_block(:, 1:nc), 1, c0, nu, nc)
2404 : END DO
2405 10 : DEALLOCATE (gram_block)
2406 :
2407 10 : CALL timestop(handle)
2408 :
2409 10 : END SUBROUTINE reduced_gram_blocks
2410 :
2411 : ! **************************************************************************************************
2412 : !> \brief fm_red_p,m+l <- red_block_pl for p <= m+nt and fm_red_m+l,p <- red_block_pl for p <= m,
2413 : !> l <= nt: the two blocks a basis extension by nt columns adds to the symmetric reduced matrix
2414 : !> \param fm_red the reduced matrix, at least m+nt rows and columns
2415 : !> \param red_block (m+nt) x nt, the new columns' overlaps with the whole extended basis
2416 : !> \param m live dimension before the extension
2417 : !> \param nt columns added
2418 : ! **************************************************************************************************
2419 60 : SUBROUTINE reduced_extend(fm_red, red_block, m, nt)
2420 :
2421 : TYPE(cp_fm_type), INTENT(IN) :: fm_red
2422 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: red_block
2423 : INTEGER, INTENT(IN) :: m, nt
2424 :
2425 60 : CALL cp_fm_set_submatrix(fm_red, red_block, 1, m + 1, m + nt, nt)
2426 60 : IF (m > 0) CALL cp_fm_set_submatrix(fm_red, red_block(1:m, :), m + 1, 1, nt, m, transpose=.TRUE.)
2427 :
2428 60 : END SUBROUTINE reduced_extend
2429 :
2430 : ! **************************************************************************************************
2431 : !> \brief A new m x m matrix on blacs_env holding the live block fm_red(1:m, 1:m)
2432 : !> \param fm_red the reduced matrix
2433 : !> \param m dimension of the live block
2434 : !> \param para_env communicator of the grid
2435 : !> \param blacs_env grid of the new matrix
2436 : !> \param name name of the new matrix
2437 : !> \param fm_live created here, released by the caller
2438 : ! **************************************************************************************************
2439 160 : SUBROUTINE reduced_live_block(fm_red, m, para_env, blacs_env, name, fm_live)
2440 :
2441 : TYPE(cp_fm_type), INTENT(IN) :: fm_red
2442 : INTEGER, INTENT(IN) :: m
2443 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2444 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2445 : CHARACTER(LEN=*), INTENT(IN) :: name
2446 : TYPE(cp_fm_type), INTENT(OUT) :: fm_live
2447 :
2448 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
2449 :
2450 80 : NULLIFY (fm_struct)
2451 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2452 80 : nrow_global=m, ncol_global=m)
2453 80 : CALL cp_fm_create(fm_live, fm_struct, name=name)
2454 80 : CALL cp_fm_struct_release(fm_struct)
2455 80 : CALL cp_fm_to_fm_submat(fm_red, fm_live, m, m, 1, 1, 1, 1)
2456 :
2457 80 : END SUBROUTINE reduced_live_block
2458 :
2459 : ! **************************************************************************************************
2460 : !> \brief fm_red_ij <- sum_pq coef_restart_pi fm_red_pq coef_restart_qj for p, q <= m and i, j <= k,
2461 : !> the rest of the matrix zeroed: two distributed products, no replicated m x m
2462 : !> \param fm_red the reduced matrix, rotated in place
2463 : !> \param m live dimension before the restart
2464 : !> \param coef_restart replicated, zero-padded to m rows
2465 : !> \param k vectors of the restart basis, the live dimension after
2466 : !> \param para_env communicator of the grid
2467 : !> \param blacs_env grid of fm_red and of the products
2468 : ! **************************************************************************************************
2469 40 : SUBROUTINE reduced_rotate(fm_red, m, coef_restart, k, para_env, blacs_env)
2470 :
2471 : TYPE(cp_fm_type), INTENT(IN) :: fm_red
2472 : INTEGER, INTENT(IN) :: m
2473 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coef_restart
2474 : INTEGER, INTENT(IN) :: k
2475 : TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2476 : TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2477 :
2478 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reduced_rotate'
2479 :
2480 : INTEGER :: handle
2481 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
2482 : TYPE(cp_fm_type) :: fm_GQ, fm_live, fm_Q, fm_QGQ
2483 :
2484 40 : CALL timeset(routineN, handle)
2485 :
2486 40 : CALL reduced_live_block(fm_red, m, para_env, blacs_env, "bse_reduced_live", fm_live)
2487 40 : NULLIFY (fm_struct)
2488 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2489 40 : nrow_global=m, ncol_global=k)
2490 40 : CALL cp_fm_create(fm_Q, fm_struct, name="bse_restart_basis")
2491 40 : CALL cp_fm_create(fm_GQ, fm_struct, name="bse_reduced_GQ")
2492 40 : CALL cp_fm_struct_release(fm_struct)
2493 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2494 40 : nrow_global=k, ncol_global=k)
2495 40 : CALL cp_fm_create(fm_QGQ, fm_struct, name="bse_reduced_QGQ")
2496 40 : CALL cp_fm_struct_release(fm_struct)
2497 :
2498 40 : CALL cp_fm_set_submatrix(fm_Q, coef_restart(1:m, 1:k))
2499 : ! GQ_pj = sum_q fm_red_pq coef_restart_qj, then QGQ_ij = sum_p coef_restart_pi GQ_pj
2500 40 : CALL parallel_gemm('N', 'N', m, k, m, 1.0_dp, fm_live, fm_Q, 0.0_dp, fm_GQ)
2501 40 : CALL parallel_gemm('T', 'N', k, k, m, 1.0_dp, fm_Q, fm_GQ, 0.0_dp, fm_QGQ)
2502 40 : CALL cp_fm_set_all(fm_red, 0.0_dp)
2503 40 : CALL cp_fm_to_fm_submat(fm_QGQ, fm_red, k, k, 1, 1, 1, 1)
2504 :
2505 40 : CALL cp_fm_release(fm_live)
2506 40 : CALL cp_fm_release(fm_Q)
2507 40 : CALL cp_fm_release(fm_GQ)
2508 40 : CALL cp_fm_release(fm_QGQ)
2509 :
2510 40 : CALL timestop(handle)
2511 :
2512 40 : END SUBROUTINE reduced_rotate
2513 :
2514 : ! **************************************************************************************************
2515 : !> \brief max_kl |sum_ia U_ia,k V_ia,l - δ_kl| over k, l <= m, in column blocks of nb
2516 : !> \param fm_U left basis
2517 : !> \param fm_V right basis, fm_U for the Euclidean case
2518 : !> \param m columns compared
2519 : !> \param nb width of one column block
2520 : !> \param para_env communicator of the row-distributed matrices
2521 : !> \return the largest deviation
2522 : ! **************************************************************************************************
2523 0 : FUNCTION orthonormality_deviation(fm_U, fm_V, m, nb, para_env) RESULT(dev)
2524 :
2525 : TYPE(cp_fm_type), INTENT(IN) :: fm_U, fm_V
2526 : INTEGER, INTENT(IN) :: m, nb
2527 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2528 : REAL(KIND=dp) :: dev
2529 :
2530 : INTEGER :: c0, l, nc
2531 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gram_block
2532 :
2533 0 : ALLOCATE (gram_block(m, nb))
2534 0 : dev = 0.0_dp
2535 0 : DO c0 = 1, m, nb
2536 0 : nc = MIN(nb, m - c0 + 1)
2537 0 : CALL subspace_gram(fm_U, 1, m, fm_V, c0, nc, para_env, gram_block(:, 1:nc))
2538 0 : DO l = 1, nc
2539 0 : gram_block(c0 + l - 1, l) = gram_block(c0 + l - 1, l) - 1.0_dp
2540 : END DO
2541 0 : dev = MAX(dev, MAXVAL(ABS(gram_block(:, 1:nc))))
2542 : END DO
2543 0 : DEALLOCATE (gram_block)
2544 :
2545 0 : END FUNCTION orthonormality_deviation
2546 :
2547 : ! **************************************************************************************************
2548 : !> \brief Euclidean norms of nt columns of a row-distributed matrix
2549 : !> Each rank sums the squares of its own rows of the npe x 1 grid, then one allreduce
2550 : !> \param fm_T row-distributed matrix
2551 : !> \param first_col first column measured
2552 : !> \param nt number of columns
2553 : !> \param para_env communicator of the rows
2554 : !> \param norms nt Euclidean norms
2555 : ! **************************************************************************************************
2556 132 : SUBROUTINE column_norms(fm_T, first_col, nt, para_env, norms)
2557 :
2558 : TYPE(cp_fm_type), INTENT(IN) :: fm_T
2559 : INTEGER, INTENT(IN) :: first_col, nt
2560 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2561 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: norms
2562 :
2563 : INTEGER :: k, nrow_local
2564 :
2565 132 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local)
2566 790 : norms(:) = 0.0_dp
2567 790 : DO k = 1, nt
2568 16582 : norms(k) = SUM(fm_T%local_data(1:nrow_local, first_col + k - 1)**2)
2569 : END DO
2570 1448 : CALL para_env%sum(norms)
2571 790 : norms(:) = SQRT(norms(:))
2572 :
2573 132 : END SUBROUTINE column_norms
2574 :
2575 : ! **************************************************************************************************
2576 : !> \brief Davidson corrections t_ia,k = r_ia,k/(d_ia - θ_k) of the selected roots k, written to fm_Z
2577 : !> Each rank divides the entries of its own rows of the npe x 1 grid, no communication
2578 : !> \param fm_R residuals, column r_offset + k belongs to root k
2579 : !> \param selected indices of the roots corrected
2580 : !> \param theta Ritz values, indexed by root
2581 : !> \param diag d_ia, replicated
2582 : !> \param fm_Z receives the corrections
2583 : !> \param first_col first column of fm_Z written
2584 : !> \param r_offset columns of fm_R before the residuals, 0 by default
2585 : ! **************************************************************************************************
2586 60 : SUBROUTINE davidson_corrections(fm_R, selected, theta, diag, fm_Z, first_col, r_offset)
2587 :
2588 : TYPE(cp_fm_type), INTENT(IN) :: fm_R
2589 : INTEGER, DIMENSION(:), INTENT(IN) :: selected
2590 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: theta, diag
2591 : TYPE(cp_fm_type), INTENT(IN) :: fm_Z
2592 : INTEGER, INTENT(IN) :: first_col
2593 : INTEGER, INTENT(IN), OPTIONAL :: r_offset
2594 :
2595 : CHARACTER(LEN=*), PARAMETER :: routineN = 'davidson_corrections'
2596 : REAL(KIND=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp)
2597 :
2598 : INTEGER :: handle, iloc, it, k, nrow_local, r_off
2599 60 : INTEGER, DIMENSION(:), POINTER :: row_indices
2600 : REAL(KIND=dp) :: denom
2601 :
2602 60 : CALL timeset(routineN, handle)
2603 :
2604 60 : r_off = 0
2605 60 : IF (PRESENT(r_offset)) r_off = r_offset
2606 60 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
2607 214 : DO it = 1, SIZE(selected)
2608 154 : k = selected(it)
2609 3910 : DO iloc = 1, nrow_local
2610 3696 : denom = diag(row_indices(iloc)) - theta(k)
2611 : ! near-degeneracy guard as in tddfpt_compute_residual_vects: threshold detects a numerically
2612 : ! zero d_ia - θ_k, which the shift (1 - eref_scale) θ_k replaces by the denominator at 0.99 θ_k
2613 3696 : IF (ABS(denom) < threshold) denom = denom + (1.0_dp - eref_scale)*theta(k)
2614 3850 : fm_Z%local_data(iloc, first_col + it - 1) = fm_R%local_data(iloc, r_off + k)/denom
2615 : END DO
2616 : END DO
2617 :
2618 60 : CALL timestop(handle)
2619 :
2620 60 : END SUBROUTINE davidson_corrections
2621 :
2622 : ! **************************************************************************************************
2623 : !> \brief T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb U_jb,m T_jb,k for m <= nv and nt columns k, twice;
2624 : !> U = V for orthonormal columns V, else the dual basis with sum_ia U_ia,m V_ia,n = δ_mn
2625 : !> One allreduce per pass inside the Gram, the subtraction stays on the rows of the npe x 1 grid
2626 : !> \param fm_T columns projected
2627 : !> \param first_col first column of fm_T projected
2628 : !> \param nt number of columns projected
2629 : !> \param fm_V the basis projected out
2630 : !> \param nv columns of fm_V
2631 : !> \param para_env communicator of the row-distributed matrices
2632 : !> \param fm_dual U
2633 : ! **************************************************************************************************
2634 64 : SUBROUTINE project_out(fm_T, first_col, nt, fm_V, nv, para_env, fm_dual)
2635 :
2636 : TYPE(cp_fm_type), INTENT(IN) :: fm_T
2637 : INTEGER, INTENT(IN) :: first_col, nt
2638 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2639 : INTEGER, INTENT(IN) :: nv
2640 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2641 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_dual
2642 :
2643 : CHARACTER(LEN=*), PARAMETER :: routineN = 'project_out'
2644 :
2645 : INTEGER :: handle, ipass, nrow_local
2646 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff
2647 :
2648 64 : IF (nt == 0 .OR. nv == 0) RETURN
2649 60 : CALL timeset(routineN, handle)
2650 60 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local)
2651 240 : ALLOCATE (coeff(nv, nt))
2652 180 : DO ipass = 1, 2
2653 120 : IF (PRESENT(fm_dual)) THEN
2654 28 : CALL subspace_gram(fm_dual, 1, nv, fm_T, first_col, nt, para_env, coeff)
2655 : ELSE
2656 92 : CALL subspace_gram(fm_V, 1, nv, fm_T, first_col, nt, para_env, coeff)
2657 : END IF
2658 180 : IF (nrow_local > 0) THEN
2659 : CALL DGEMM('N', 'N', nrow_local, nt, nv, -1.0_dp, &
2660 : fm_V%local_data, SIZE(fm_V%local_data, 1), coeff, nv, 1.0_dp, &
2661 120 : fm_T%local_data(:, first_col:first_col + nt - 1), SIZE(fm_T%local_data, 1))
2662 : END IF
2663 : END DO
2664 60 : DEALLOCATE (coeff)
2665 60 : CALL timestop(handle)
2666 :
2667 60 : END SUBROUTINE project_out
2668 :
2669 : ! **************************************************************************************************
2670 : !> \brief Normalises nt columns and removes those that vanished in the projection
2671 : !> Only the norms communicate; the scaling and the zero-fill stay on the local rows
2672 : !> \param fm_T row-distributed matrix; the kept columns move to the front of the block
2673 : !> \param first_col first column of the block
2674 : !> \param nt on exit the number of columns kept
2675 : !> \param para_env communicator of the rows
2676 : ! **************************************************************************************************
2677 64 : SUBROUTINE drop_small_columns(fm_T, first_col, nt, para_env)
2678 :
2679 : TYPE(cp_fm_type), INTENT(IN) :: fm_T
2680 : INTEGER, INTENT(IN) :: first_col
2681 : INTEGER, INTENT(INOUT) :: nt
2682 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2683 :
2684 : REAL(KIND=dp), PARAMETER :: col_norm_drop = 1.0E-10_dp
2685 :
2686 : INTEGER :: k, n_kept, nrow_local
2687 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: norms
2688 :
2689 64 : IF (nt == 0) RETURN
2690 64 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local)
2691 192 : ALLOCATE (norms(nt))
2692 64 : CALL column_norms(fm_T, first_col, nt, para_env, norms)
2693 64 : n_kept = 0
2694 246 : DO k = 1, nt
2695 : ! absolute: the corrections r/(d - θ) arrive unnormalised, a column below this is numerically zero
2696 182 : IF (norms(k) < col_norm_drop) CYCLE
2697 182 : n_kept = n_kept + 1
2698 : fm_T%local_data(1:nrow_local, first_col + n_kept - 1) = &
2699 4614 : fm_T%local_data(1:nrow_local, first_col + k - 1)/norms(k)
2700 : END DO
2701 64 : fm_T%local_data(1:nrow_local, first_col + n_kept:first_col + nt - 1) = 0.0_dp
2702 64 : nt = n_kept
2703 64 : DEALLOCATE (norms)
2704 :
2705 64 : END SUBROUTINE drop_small_columns
2706 :
2707 : ! **************************************************************************************************
2708 : !> \brief Orthonormalises nt columns by two Cholesky QR passes on the replicated Gram matrix,
2709 : !> T_ia,k <- sum_l T_ia,l (R^-1)_lk with sum_p R_pk R_pl = G_kl = sum_ia T_ia,k T_ia,l. A
2710 : !> failed factorization means a numerically dependent block: the directions of the Gram matrix
2711 : !> below the drop tolerance are removed (canonical orthogonalisation), nt shrinks, and the
2712 : !> passes start again. The Gram is the only communication: R is factorised on one rank and
2713 : !> broadcast, the triangular solve runs on the rows of the npe x 1 grid
2714 : !> \param fm_T the block to orthonormalise, columns first_col .. first_col+nt-1
2715 : !> \param first_col first column of the block
2716 : !> \param nt columns of the block; on exit the number kept
2717 : !> \param para_env communicator of the row-distributed matrix
2718 : !> \param ok .FALSE. when two drop rounds still leave no factorisable Gram matrix, or no direction survives
2719 : !> \param n_dropped columns removed as dependent
2720 : ! **************************************************************************************************
2721 64 : SUBROUTINE cholqr2(fm_T, first_col, nt, para_env, ok, n_dropped)
2722 :
2723 : TYPE(cp_fm_type), INTENT(IN) :: fm_T
2724 : INTEGER, INTENT(IN) :: first_col
2725 : INTEGER, INTENT(INOUT) :: nt
2726 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2727 : LOGICAL, INTENT(OUT) :: ok
2728 : INTEGER, INTENT(OUT), OPTIONAL :: n_dropped
2729 :
2730 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholqr2'
2731 : REAL(KIND=dp), PARAMETER :: eig_drop_rel = 1.0E-10_dp
2732 :
2733 : INTEGER :: handle, info, ipass, k, n_drop_calls, &
2734 : n_keep, nrow_local
2735 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig
2736 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chol_factor, gram, gram_eigvec, T_reduced
2737 :
2738 64 : ok = .TRUE.
2739 64 : IF (PRESENT(n_dropped)) n_dropped = 0
2740 64 : IF (nt == 0) RETURN
2741 64 : CALL timeset(routineN, handle)
2742 64 : CALL cp_fm_get_info(fm_T, nrow_local=nrow_local)
2743 384 : ALLOCATE (gram(nt, nt), chol_factor(nt, nt))
2744 :
2745 64 : n_drop_calls = 0
2746 64 : ipass = 0
2747 192 : DO WHILE (ipass < 2)
2748 128 : ipass = ipass + 1
2749 : ! G_kl = sum_ia T_ia,k T_ia,l = sum_p R_pk R_pl
2750 128 : CALL subspace_gram(fm_T, first_col, nt, fm_T, first_col, nt, para_env, gram)
2751 1840 : chol_factor(:, :) = gram(:, :)
2752 128 : info = 0
2753 128 : IF (para_env%is_source()) CALL DPOTRF('U', nt, chol_factor, nt, info)
2754 128 : CALL para_env%bcast(info)
2755 128 : IF (info /= 0) THEN
2756 0 : n_drop_calls = n_drop_calls + 1
2757 0 : IF (n_drop_calls > 2) THEN
2758 0 : ok = .FALSE.
2759 0 : EXIT
2760 : END IF
2761 : ! T_ia,k <- sum_l T_ia,l U_lk Λ_k^-1/2 over the eigenpairs of G above the drop tolerance,
2762 : ! eig_drop_rel of the largest eigenvalue: directions below it are numerically dependent
2763 0 : ALLOCATE (eig(nt), gram_eigvec(nt, nt))
2764 0 : CALL solve_replicated(gram, nt, para_env, eig, gram_eigvec)
2765 0 : n_keep = COUNT(eig > eig_drop_rel*eig(nt))
2766 0 : IF (n_keep == 0) THEN
2767 0 : ok = .FALSE.
2768 0 : DEALLOCATE (eig, gram_eigvec)
2769 0 : EXIT
2770 : END IF
2771 0 : DO k = 1, n_keep
2772 0 : gram_eigvec(:, nt - n_keep + k) = gram_eigvec(:, nt - n_keep + k)/SQRT(eig(nt - n_keep + k))
2773 : END DO
2774 0 : IF (nrow_local > 0) THEN
2775 0 : ALLOCATE (T_reduced(nrow_local, n_keep))
2776 : CALL DGEMM('N', 'N', nrow_local, n_keep, nt, 1.0_dp, &
2777 : fm_T%local_data(:, first_col:first_col + nt - 1), SIZE(fm_T%local_data, 1), &
2778 0 : gram_eigvec(:, nt - n_keep + 1:nt), nt, 0.0_dp, T_reduced, nrow_local)
2779 0 : fm_T%local_data(1:nrow_local, first_col:first_col + n_keep - 1) = T_reduced(:, :)
2780 0 : fm_T%local_data(1:nrow_local, first_col + n_keep:first_col + nt - 1) = 0.0_dp
2781 0 : DEALLOCATE (T_reduced)
2782 : END IF
2783 0 : IF (PRESENT(n_dropped)) n_dropped = n_dropped + nt - n_keep
2784 0 : nt = n_keep
2785 0 : DEALLOCATE (eig, gram_eigvec, gram, chol_factor)
2786 0 : ALLOCATE (gram(nt, nt), chol_factor(nt, nt))
2787 0 : ipass = 0
2788 0 : CYCLE
2789 : END IF
2790 : ! T_ia,k <- sum_l T_ia,l (R^-1)_lk on the local rows
2791 128 : CALL para_env%bcast(chol_factor)
2792 192 : IF (nrow_local > 0) THEN
2793 : CALL DTRSM('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
2794 128 : fm_T%local_data(:, first_col:first_col + nt - 1), SIZE(fm_T%local_data, 1))
2795 : END IF
2796 : END DO
2797 :
2798 64 : DEALLOCATE (gram, chol_factor)
2799 64 : CALL timestop(handle)
2800 :
2801 128 : END SUBROUTINE cholqr2
2802 :
2803 : ! **************************************************************************************************
2804 : !> \brief Makes the nt columns T after column m of fm_V orthonormal to the basis and to each other:
2805 : !> T_ia,k <- T_ia,k - sum_p V_ia,p sum_jb U_jb,p T_jb,k twice, the columns that vanished dropped, then
2806 : !> the Euclidean cholqr2; nt shrinks by the vanished and the dependent columns
2807 : !> \param fm_V basis in columns 1..m, the block T in columns m+1..m+nt
2808 : !> \param m columns of the basis
2809 : !> \param nt columns of T; on exit the number kept
2810 : !> \param para_env communicator of the row-distributed matrix
2811 : !> \param n_dependent counter of the columns dropped as dependent, incremented
2812 : !> \param fm_dual U, the dual basis with sum_ia U_ia,p V_ia,q = δ_pq; fm_V itself when absent
2813 : ! **************************************************************************************************
2814 64 : SUBROUTINE extend_orthonormal(fm_V, m, nt, para_env, n_dependent, fm_dual)
2815 :
2816 : TYPE(cp_fm_type), INTENT(IN) :: fm_V
2817 : INTEGER, INTENT(IN) :: m
2818 : INTEGER, INTENT(INOUT) :: nt
2819 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2820 : INTEGER, INTENT(INOUT) :: n_dependent
2821 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_dual
2822 :
2823 : INTEGER :: n_dropped
2824 : LOGICAL :: ok
2825 :
2826 64 : CALL project_out(fm_V, m + 1, nt, fm_V, m, para_env, fm_dual)
2827 64 : CALL drop_small_columns(fm_V, m + 1, nt, para_env)
2828 64 : CALL cholqr2(fm_V, m + 1, nt, para_env, ok, n_dropped)
2829 : ! cholqr2 gives up after two rounds of dropping dependent directions, or when none is left
2830 64 : IF (.NOT. ok) CPABORT("BSE Davidson: orthonormalisation broke down")
2831 64 : n_dependent = n_dependent + n_dropped
2832 :
2833 64 : END SUBROUTINE extend_orthonormal
2834 :
2835 : ! **************************************************************************************************
2836 : !> \brief Largest deviation of the nt new columns from orthonormality against the whole basis,
2837 : !> sum_ia V_ia,m D_ia,n - δ_mn for m = 1..nv+nt, n = nv+1..nv+nt; D = V in the Euclidean case
2838 : !> \param fm_V the basis with its nt new columns
2839 : !> \param fm_dual D, fm_V itself in the Euclidean case
2840 : !> \param nv columns of the basis before the extension
2841 : !> \param nt new columns
2842 : !> \param para_env communicator of the row-distributed matrices
2843 : !> \return the largest deviation
2844 : ! **************************************************************************************************
2845 66 : FUNCTION extension_deviation(fm_V, fm_dual, nv, nt, para_env) RESULT(dev)
2846 :
2847 : TYPE(cp_fm_type), INTENT(IN) :: fm_V, fm_dual
2848 : INTEGER, INTENT(IN) :: nv, nt
2849 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2850 : REAL(KIND=dp) :: dev
2851 :
2852 : INTEGER :: k
2853 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gram
2854 :
2855 66 : dev = 0.0_dp
2856 66 : IF (nt == 0) RETURN
2857 264 : ALLOCATE (gram(nv + nt, nt))
2858 66 : CALL subspace_gram(fm_V, 1, nv + nt, fm_dual, nv + 1, nt, para_env, gram)
2859 262 : DO k = 1, nt
2860 262 : gram(nv + k, k) = gram(nv + k, k) - 1.0_dp
2861 : END DO
2862 3068 : dev = MAXVAL(ABS(gram))
2863 66 : DEALLOCATE (gram)
2864 :
2865 66 : END FUNCTION extension_deviation
2866 :
2867 : ! **************************************************************************************************
2868 : !> \brief Debug check of a Davidson result against the full diagonalization of the explicit
2869 : !> matrices: energies, and per multiplet the smallest singular value of X_ref^T X for the
2870 : !> TDA, or of X_ref^T X - Y_ref^T Y for ABBA. ABBA diagonalizes the Hermitian reduction
2871 : !> C = (A-B)^1/2 (A+B) (A-B)^1/2 with C T^n = (Ω^n)^2 T^n,
2872 : !> (X+Y)^n = (Ω^n)^-1/2 (A-B)^1/2 T^n, (X-Y)^n = (Ω^n)^1/2 (A-B)^-1/2 T^n.
2873 : !> A failed reference solve aborts: this is debug output with nothing to fall back on
2874 : !> \param fm_A_explicit A from create_A_and_B, N_ov x N_ov
2875 : !> \param exc_ens Davidson energies in Hartree
2876 : !> \param fm_X Davidson X_ia^n, column n
2877 : !> \param mp2_env handed to create_hermitian_form_of_ABBA
2878 : !> \param unit_nr output unit, positive on the writing rank only
2879 : !> \param fm_B_explicit present for an ABBA result, together with fm_Y
2880 : !> \param fm_Y Davidson Y_ia^n, with fm_B_explicit
2881 : ! **************************************************************************************************
2882 0 : SUBROUTINE bse_davidson_refcheck(fm_A_explicit, exc_ens, fm_X, mp2_env, unit_nr, fm_B_explicit, fm_Y)
2883 :
2884 : TYPE(cp_fm_type), INTENT(IN) :: fm_A_explicit
2885 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_ens
2886 : TYPE(cp_fm_type), INTENT(IN) :: fm_X
2887 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2888 : INTEGER, INTENT(IN) :: unit_nr
2889 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_B_explicit, fm_Y
2890 :
2891 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_davidson_refcheck'
2892 :
2893 : INTEGER :: diag_info, handle, k, n_ov, n_ref, n_want
2894 : REAL(KIND=dp) :: min_overlap
2895 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ref_ens
2896 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: diff_ref, sum_ref, X_dav, X_ref, Y_dav, &
2897 0 : Y_ref
2898 : TYPE(cp_fm_type) :: fm_C, fm_eigvec, fm_inv_sqrt_A_minus_B, &
2899 : fm_sqrt_A_minus_B, fm_work
2900 :
2901 0 : CALL timeset(routineN, handle)
2902 :
2903 : ! B and Y arrive together or not at all
2904 0 : CPASSERT(PRESENT(fm_B_explicit) .EQV. PRESENT(fm_Y))
2905 0 : IF (unit_nr > 0) THEN
2906 0 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
2907 0 : IF (PRESENT(fm_B_explicit)) THEN
2908 0 : WRITE (unit_nr, '(T2,A10,T13,A)') 'BSE|DEBUG|', 'Reference check against the explicit A and B'
2909 : ELSE
2910 0 : WRITE (unit_nr, '(T2,A10,T13,A)') 'BSE|DEBUG|', 'Reference check against the explicit A'
2911 : END IF
2912 : END IF
2913 0 : CALL cp_fm_get_info(fm_A_explicit, nrow_global=n_ov)
2914 0 : n_want = SIZE(exc_ens)
2915 0 : ALLOCATE (ref_ens(n_ov))
2916 :
2917 0 : IF (PRESENT(fm_B_explicit)) THEN
2918 : CALL create_hermitian_form_of_ABBA(fm_A_explicit, fm_B_explicit, fm_C, fm_sqrt_A_minus_B, &
2919 0 : fm_inv_sqrt_A_minus_B, unit_nr, mp2_env, 0.0_dp)
2920 0 : CALL cp_fm_create(fm_eigvec, fm_C%matrix_struct)
2921 0 : CALL choose_eigv_solver(fm_C, fm_eigvec, ref_ens, diag_info)
2922 0 : IF (diag_info /= 0) CPABORT("Reference diagonalization of C failed in the BSE Davidson check")
2923 0 : IF (ref_ens(1) <= 0.0_dp) THEN
2924 : CALL cp_abort(__LOCATION__, &
2925 0 : "Reference matrix C has a non-positive eigenvalue in the BSE Davidson check")
2926 : END IF
2927 0 : ref_ens(:) = SQRT(ref_ens(:))
2928 :
2929 0 : n_ref = multiplet_end(ref_ens, n_want, n_ov, deg_thresh)
2930 0 : ALLOCATE (sum_ref(n_ov, n_ref), diff_ref(n_ov, n_ref))
2931 0 : CALL cp_fm_create(fm_work, fm_C%matrix_struct)
2932 : CALL parallel_gemm("N", "N", n_ov, n_ref, n_ov, 1.0_dp, fm_sqrt_A_minus_B, fm_eigvec, 0.0_dp, &
2933 0 : fm_work)
2934 0 : CALL cp_fm_get_submatrix(fm_work, sum_ref, 1, 1, n_ov, n_ref)
2935 : CALL parallel_gemm("N", "N", n_ov, n_ref, n_ov, 1.0_dp, fm_inv_sqrt_A_minus_B, fm_eigvec, &
2936 0 : 0.0_dp, fm_work)
2937 0 : CALL cp_fm_get_submatrix(fm_work, diff_ref, 1, 1, n_ov, n_ref)
2938 0 : CALL cp_fm_release(fm_work)
2939 0 : CALL cp_fm_release(fm_C)
2940 0 : CALL cp_fm_release(fm_sqrt_A_minus_B)
2941 0 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
2942 :
2943 0 : ALLOCATE (X_ref(n_ov, n_ref), Y_ref(n_ov, n_ref), Y_dav(n_ov, n_want))
2944 0 : DO k = 1, n_ref
2945 0 : sum_ref(:, k) = sum_ref(:, k)/SQRT(ref_ens(k))
2946 0 : diff_ref(:, k) = diff_ref(:, k)*SQRT(ref_ens(k))
2947 : END DO
2948 0 : X_ref(:, :) = 0.5_dp*(sum_ref(:, :) + diff_ref(:, :))
2949 0 : Y_ref(:, :) = 0.5_dp*(sum_ref(:, :) - diff_ref(:, :))
2950 0 : DEALLOCATE (sum_ref, diff_ref)
2951 0 : CALL cp_fm_get_submatrix(fm_Y, Y_dav)
2952 : ELSE
2953 : ! the eigensolver destroys its input
2954 0 : CALL cp_fm_create(fm_work, fm_A_explicit%matrix_struct)
2955 0 : CALL cp_fm_to_fm(fm_A_explicit, fm_work)
2956 0 : CALL cp_fm_create(fm_eigvec, fm_A_explicit%matrix_struct)
2957 0 : CALL choose_eigv_solver(fm_work, fm_eigvec, ref_ens, diag_info)
2958 0 : IF (diag_info /= 0) CPABORT("Reference diagonalization of A failed in the BSE Davidson check")
2959 0 : CALL cp_fm_release(fm_work)
2960 :
2961 0 : n_ref = multiplet_end(ref_ens, n_want, n_ov, deg_thresh)
2962 0 : ALLOCATE (X_ref(n_ov, n_ref))
2963 0 : CALL cp_fm_get_submatrix(fm_eigvec, X_ref, 1, 1, n_ov, n_ref)
2964 : END IF
2965 0 : CALL cp_fm_release(fm_eigvec)
2966 :
2967 0 : ALLOCATE (X_dav(n_ov, n_want))
2968 0 : CALL cp_fm_get_submatrix(fm_X, X_dav)
2969 :
2970 : ! Y_ref and Y_dav stay unallocated for the TDA and then count as absent (Fortran 2008)
2971 0 : CALL min_multiplet_overlap(ref_ens, n_want, n_ref, X_ref, X_dav, min_overlap, Y_ref, Y_dav)
2972 0 : CALL print_refcheck(MAXVAL(ABS(exc_ens(:) - ref_ens(1:n_want))), min_overlap, unit_nr)
2973 :
2974 0 : DEALLOCATE (ref_ens, X_ref, X_dav)
2975 0 : IF (ALLOCATED(Y_ref)) DEALLOCATE (Y_ref, Y_dav)
2976 :
2977 0 : CALL timestop(handle)
2978 :
2979 0 : END SUBROUTINE bse_davidson_refcheck
2980 :
2981 : ! **************************************************************************************************
2982 : !> \brief Smallest singular value over the multiplets of O_kl = sum_ia (X_ref,ia^k X_ia^l - Y_ref,ia^k Y_ia^l),
2983 : !> k over a multiplet of the reference and l over the Davidson states inside it. Degenerate
2984 : !> eigenvectors are arbitrary within their multiplet, so the vectors are not compared one by one:
2985 : !> every singular value of O is 1 when the Davidson states lie in the span of the multiplet
2986 : !> \param ref_ens reference energies in Hartree, ascending
2987 : !> \param n_want Davidson states compared
2988 : !> \param n_ref reference states up to the end of the multiplet that contains state n_want
2989 : !> \param X_ref reference X, N_ov x n_ref
2990 : !> \param X_dav Davidson X, N_ov x n_want
2991 : !> \param min_overlap the smallest singular value found, 1 if every multiplet is reproduced
2992 : !> \param Y_ref reference Y for ABBA, with Y_dav
2993 : !> \param Y_dav Davidson Y for ABBA
2994 : ! **************************************************************************************************
2995 0 : SUBROUTINE min_multiplet_overlap(ref_ens, n_want, n_ref, X_ref, X_dav, min_overlap, Y_ref, Y_dav)
2996 :
2997 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: ref_ens
2998 : INTEGER, INTENT(IN) :: n_want, n_ref
2999 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: X_ref, X_dav
3000 : REAL(KIND=dp), INTENT(OUT) :: min_overlap
3001 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
3002 : OPTIONAL :: Y_ref, Y_dav
3003 :
3004 : INTEGER :: info, lwork, mult_first, mult_last, &
3005 : mult_last_dav
3006 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_vals, work
3007 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: overlap
3008 : REAL(KIND=dp), DIMENSION(1, 1) :: dummy
3009 :
3010 0 : min_overlap = 1.0_dp
3011 0 : mult_first = 1
3012 0 : DO WHILE (mult_first <= n_want)
3013 0 : mult_last = multiplet_end(ref_ens, mult_first, n_ref, deg_thresh)
3014 0 : mult_last_dav = MIN(mult_last, n_want)
3015 0 : ALLOCATE (overlap(mult_last - mult_first + 1, mult_last_dav - mult_first + 1), sing_vals(mult_last_dav - mult_first + 1))
3016 0 : overlap(:, :) = MATMUL(TRANSPOSE(X_ref(:, mult_first:mult_last)), X_dav(:, mult_first:mult_last_dav))
3017 0 : IF (PRESENT(Y_ref)) THEN
3018 0 : overlap(:, :) = overlap(:, :) - MATMUL(TRANSPOSE(Y_ref(:, mult_first:mult_last)), Y_dav(:, mult_first:mult_last_dav))
3019 : END IF
3020 : ! DGESVD wants lwork >= MAX(3 n + m, 5 n) for the m x n overlap with n <= m; 5 m + 10 covers both
3021 0 : lwork = 5*(mult_last - mult_first + 1) + 10
3022 0 : ALLOCATE (work(lwork))
3023 : CALL DGESVD('N', 'N', mult_last - mult_first + 1, mult_last_dav - mult_first + 1, overlap, mult_last - mult_first + 1, sing_vals, &
3024 0 : dummy, 1, dummy, 1, work, lwork, info)
3025 0 : IF (info /= 0) CPABORT("SVD failed in the BSE Davidson check")
3026 0 : min_overlap = MIN(min_overlap, MINVAL(sing_vals))
3027 0 : DEALLOCATE (overlap, sing_vals, work)
3028 0 : mult_first = mult_last + 1
3029 : END DO
3030 :
3031 0 : END SUBROUTINE min_multiplet_overlap
3032 :
3033 : ! **************************************************************************************************
3034 : !> \brief The two lines of a reference check
3035 : !> \param dev_E largest energy deviation in Hartree
3036 : !> \param min_overlap smallest singular value from min_multiplet_overlap
3037 : !> \param unit_nr output unit, positive on the writing rank only
3038 : ! **************************************************************************************************
3039 0 : SUBROUTINE print_refcheck(dev_E, min_overlap, unit_nr)
3040 :
3041 : REAL(KIND=dp), INTENT(IN) :: dev_E, min_overlap
3042 : INTEGER, INTENT(IN) :: unit_nr
3043 :
3044 0 : IF (unit_nr > 0) THEN
3045 0 : WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
3046 0 : 'Max energy deviation vs full diagonalization (eV)', dev_E*evolt
3047 0 : WRITE (unit_nr, '(T2,A10,T13,A,T63,F18.12)') 'BSE|DEBUG|', &
3048 0 : 'Min subspace overlap vs full diagonalization', min_overlap
3049 : END IF
3050 :
3051 0 : END SUBROUTINE print_refcheck
3052 :
3053 0 : END MODULE bse_davidson
|
