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 module that contains the algorithms to perform an iterative
10 : !> diagonalization by the block-Davidson approach
11 : !> P. Blaha, et al J. Comp. Physics, 229, (2010), 453-460
12 : !> Iterative diagonalization in augmented plane wave based
13 : !> methods in electronic structure calculations
14 : !> \par History
15 : !> 05.2011 created [MI]
16 : !> \author MI
17 : ! **************************************************************************************************
18 : MODULE qs_scf_block_davidson
19 :
20 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
21 : cp_cfm_gemm,&
22 : cp_cfm_get_diag,&
23 : cp_cfm_scale_and_add,&
24 : cp_cfm_vectorsnorm
25 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_get_submatrix,&
28 : cp_cfm_release,&
29 : cp_cfm_set_all,&
30 : cp_cfm_set_submatrix,&
31 : cp_cfm_to_cfm,&
32 : cp_cfm_type
33 : USE cp_dbcsr_api, ONLY: &
34 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_init_p, &
35 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
36 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release_p, dbcsr_type, &
37 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
38 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_diag,&
39 : dbcsr_scale_by_vector
40 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 : copy_fm_to_dbcsr,&
42 : cp_dbcsr_m_by_n_from_row_template,&
43 : cp_dbcsr_m_by_n_from_template,&
44 : cp_dbcsr_sm_fm_multiply
45 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
46 : cp_fm_scale_and_add,&
47 : cp_fm_symm,&
48 : cp_fm_transpose,&
49 : cp_fm_triangular_invert,&
50 : cp_fm_uplo_to_full
51 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
52 : cp_fm_cholesky_restore
53 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
54 : cp_fm_power
55 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
56 : cp_fm_struct_release,&
57 : cp_fm_struct_type
58 : USE cp_fm_types, ONLY: cp_fm_create,&
59 : cp_fm_get_diag,&
60 : cp_fm_release,&
61 : cp_fm_set_all,&
62 : cp_fm_to_fm,&
63 : cp_fm_to_fm_submat,&
64 : cp_fm_type,&
65 : cp_fm_vectorsnorm
66 : USE kinds, ONLY: dp
67 : USE machine, ONLY: m_walltime
68 : USE mathconstants, ONLY: z_one,&
69 : z_zero
70 : USE message_passing, ONLY: mp_comm_type
71 : USE parallel_gemm_api, ONLY: parallel_gemm
72 : USE preconditioner, ONLY: apply_preconditioner
73 : USE preconditioner_types, ONLY: preconditioner_type
74 : USE qs_block_davidson_types, ONLY: davidson_type
75 : USE qs_mo_types, ONLY: get_mo_set,&
76 : mo_set_type
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 : PRIVATE
81 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_block_davidson'
82 :
83 : PUBLIC :: generate_extended_space, generate_extended_space_c, generate_extended_space_sparse
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param bdav_env ...
90 : !> \param mo_set ...
91 : !> \param matrix_h ...
92 : !> \param matrix_s ...
93 : !> \param output_unit ...
94 : !> \param preconditioner ...
95 : ! **************************************************************************************************
96 40 : SUBROUTINE generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
97 : preconditioner)
98 :
99 : TYPE(davidson_type) :: bdav_env
100 : TYPE(mo_set_type), INTENT(IN) :: mo_set
101 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
102 : INTEGER, INTENT(IN) :: output_unit
103 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
104 :
105 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space'
106 :
107 : INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, max_iter, n, nao, nmat, nmat2, &
108 : nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
109 40 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
110 40 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
111 : LOGICAL :: converged, do_apply_preconditioner
112 : REAL(dp) :: lambda, max_norm, min_norm, t1, t2
113 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm
114 40 : REAL(dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals
115 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
116 : TYPE(cp_fm_type) :: c_conv, c_notconv, c_out, h_block, h_fm, &
117 : m_hc, m_sc, m_tmp, mt_tmp, s_block, &
118 : s_fm, v_block, w_block
119 : TYPE(cp_fm_type), POINTER :: c_pz, c_z, mo_coeff
120 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
121 :
122 40 : CALL timeset(routineN, handle)
123 :
124 40 : NULLIFY (mo_coeff, mo_coeff_b, eigenvalues)
125 :
126 40 : do_apply_preconditioner = .FALSE.
127 40 : IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE.
128 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, eigenvalues=eigenvalues, &
129 40 : nao=nao, nmo=nmo, homo=homo)
130 40 : IF (do_apply_preconditioner) THEN
131 36 : max_iter = bdav_env%max_iter
132 : ELSE
133 : max_iter = 1
134 : END IF
135 :
136 40 : NULLIFY (c_z, c_pz)
137 40 : NULLIFY (evals, eig_not_conv)
138 40 : t1 = m_walltime()
139 40 : IF (output_unit > 0) THEN
140 : WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
141 0 : " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
142 : END IF
143 :
144 120 : ALLOCATE (iconv(nmo))
145 80 : ALLOCATE (inotconv(nmo))
146 120 : ALLOCATE (ritz_coeff(nmo))
147 80 : ALLOCATE (vnorm(nmo))
148 :
149 40 : converged = .FALSE.
150 124 : DO iter = 1, max_iter
151 :
152 : ! compute Ritz values
153 88 : ritz_coeff = 0.0_dp
154 88 : CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc")
155 88 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo)
156 88 : CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc")
157 88 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_coeff, m_sc, nmo)
158 :
159 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
160 : context=mo_coeff%matrix_struct%context, &
161 88 : para_env=mo_coeff%matrix_struct%para_env)
162 88 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
163 88 : CALL cp_fm_struct_release(fm_struct_tmp)
164 :
165 88 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp)
166 88 : CALL cp_fm_get_diag(m_tmp, ritz_coeff)
167 88 : CALL cp_fm_release(m_tmp)
168 :
169 : ! Check for converged eigenvectors
170 88 : c_z => bdav_env%matrix_z
171 88 : c_pz => bdav_env%matrix_pz
172 88 : CALL cp_fm_to_fm(m_sc, c_z)
173 88 : CALL cp_fm_column_scale(c_z, ritz_coeff)
174 88 : CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc)
175 88 : CALL cp_fm_vectorsnorm(c_z, vnorm)
176 :
177 88 : nmo_converged = 0
178 88 : nmo_not_converged = 0
179 88 : max_norm = 0.0_dp
180 88 : min_norm = 1.e10_dp
181 3256 : DO imo = 1, nmo
182 3168 : max_norm = MAX(max_norm, vnorm(imo))
183 3256 : min_norm = MIN(min_norm, vnorm(imo))
184 : END DO
185 88 : iconv = 0
186 88 : inotconv = 0
187 3256 : DO imo = 1, nmo
188 3256 : IF (vnorm(imo) <= bdav_env%eps_iter) THEN
189 136 : nmo_converged = nmo_converged + 1
190 136 : iconv(nmo_converged) = imo
191 : ELSE
192 3032 : nmo_not_converged = nmo_not_converged + 1
193 3032 : inotconv(nmo_not_converged) = imo
194 : END IF
195 : END DO
196 :
197 88 : IF (nmo_converged > 0) THEN
198 48 : ALLOCATE (iconv_set(nmo_converged, 2))
199 48 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
200 16 : i_last = iconv(1)
201 16 : nset = 0
202 152 : DO j = 1, nmo_converged
203 136 : imo = iconv(j)
204 :
205 152 : IF (imo == i_last + 1) THEN
206 102 : i_last = imo
207 102 : iconv_set(nset, 2) = imo
208 : ELSE
209 34 : i_last = imo
210 34 : nset = nset + 1
211 34 : iconv_set(nset, 1) = imo
212 34 : iconv_set(nset, 2) = imo
213 : END IF
214 : END DO
215 16 : nset_conv = nset
216 :
217 16 : i_last = inotconv(1)
218 16 : nset = 0
219 456 : DO j = 1, nmo_not_converged
220 440 : imo = inotconv(j)
221 :
222 456 : IF (imo == i_last + 1) THEN
223 398 : i_last = imo
224 398 : inotconv_set(nset, 2) = imo
225 : ELSE
226 42 : i_last = imo
227 42 : nset = nset + 1
228 42 : inotconv_set(nset, 1) = imo
229 42 : inotconv_set(nset, 2) = imo
230 : END IF
231 : END DO
232 16 : nset_not_conv = nset
233 16 : CALL cp_fm_release(m_sc)
234 16 : CALL cp_fm_release(m_hc)
235 16 : NULLIFY (c_z, c_pz)
236 : END IF
237 :
238 88 : IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
239 4 : converged = .TRUE.
240 4 : DEALLOCATE (iconv_set)
241 4 : DEALLOCATE (inotconv_set)
242 4 : t2 = m_walltime()
243 4 : IF (output_unit > 0) THEN
244 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
245 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
246 :
247 0 : WRITE (output_unit, *) " Reached convergence in ", iter, &
248 0 : " Davidson iterations"
249 : END IF
250 :
251 : EXIT
252 : END IF
253 :
254 84 : IF (nmo_converged > 0) THEN
255 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
256 : context=mo_coeff%matrix_struct%context, &
257 12 : para_env=mo_coeff%matrix_struct%para_env)
258 : !allocate h_fm
259 12 : CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp")
260 : !allocate s_fm
261 12 : CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp")
262 : !copy matrix_h in h_fm
263 12 : CALL copy_dbcsr_to_fm(matrix_h, h_fm)
264 12 : CALL cp_fm_uplo_to_full(h_fm, s_fm)
265 :
266 : !copy matrix_s in s_fm
267 : ! CALL cp_fm_set_all(s_fm,0.0_dp)
268 12 : CALL copy_dbcsr_to_fm(matrix_s, s_fm)
269 :
270 : !allocate c_out
271 12 : CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp")
272 : ! set c_out to zero
273 12 : CALL cp_fm_set_all(c_out, 0.0_dp)
274 12 : CALL cp_fm_struct_release(fm_struct_tmp)
275 :
276 : !allocate c_conv
277 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
278 : context=mo_coeff%matrix_struct%context, &
279 12 : para_env=mo_coeff%matrix_struct%para_env)
280 12 : CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv")
281 12 : CALL cp_fm_set_all(c_conv, 0.0_dp)
282 : !allocate m_tmp
283 12 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc")
284 12 : CALL cp_fm_struct_release(fm_struct_tmp)
285 : END IF
286 :
287 : !allocate c_notconv
288 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
289 : context=mo_coeff%matrix_struct%context, &
290 84 : para_env=mo_coeff%matrix_struct%para_env)
291 84 : CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv")
292 84 : CALL cp_fm_set_all(c_notconv, 0.0_dp)
293 84 : IF (nmo_converged > 0) THEN
294 12 : CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc")
295 12 : CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc")
296 : !allocate c_z
297 12 : ALLOCATE (c_z, c_pz)
298 12 : CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z")
299 12 : CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz")
300 12 : CALL cp_fm_set_all(c_z, 0.0_dp)
301 :
302 : ! sum contributions to c_out
303 12 : jj = 1
304 34 : DO j = 1, nset_conv
305 22 : i_first = iconv_set(j, 1)
306 22 : i_last = iconv_set(j, 2)
307 22 : n = i_last - i_first + 1
308 22 : CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj)
309 34 : jj = jj + n
310 : END DO
311 12 : CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp)
312 12 : CALL parallel_gemm('N', 'T', nao, nao, nmo_converged, 1.0_dp, m_tmp, m_tmp, 0.0_dp, c_out)
313 :
314 : ! project c_out out of H
315 12 : lambda = 100.0_dp*ABS(eigenvalues(homo))
316 12 : CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm)
317 12 : CALL cp_fm_release(m_tmp)
318 12 : CALL cp_fm_release(h_fm)
319 :
320 : END IF
321 :
322 : !allocate m_tmp
323 84 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm")
324 84 : CALL cp_fm_struct_release(fm_struct_tmp)
325 84 : IF (nmo_converged > 0) THEN
326 36 : ALLOCATE (eig_not_conv(nmo_not_converged))
327 12 : jj = 1
328 42 : DO j = 1, nset_not_conv
329 30 : i_first = inotconv_set(j, 1)
330 30 : i_last = inotconv_set(j, 2)
331 30 : n = i_last - i_first + 1
332 30 : CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj)
333 404 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
334 42 : jj = jj + n
335 : END DO
336 12 : CALL parallel_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc)
337 12 : CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_notconv, 0.0_dp, m_sc)
338 : ! extend suspace using only the not converged vectors
339 12 : CALL cp_fm_to_fm(m_sc, m_tmp)
340 12 : CALL cp_fm_column_scale(m_tmp, eig_not_conv)
341 12 : CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc)
342 12 : DEALLOCATE (eig_not_conv)
343 12 : CALL cp_fm_to_fm(m_tmp, c_z)
344 : ELSE
345 72 : CALL cp_fm_to_fm(mo_coeff, c_notconv)
346 : END IF
347 :
348 : !preconditioner
349 84 : IF (do_apply_preconditioner) THEN
350 80 : IF (preconditioner%in_use /= 0) THEN
351 80 : CALL apply_preconditioner(preconditioner, c_z, c_pz)
352 : ELSE
353 0 : CALL cp_fm_to_fm(c_z, c_pz)
354 : END IF
355 : ELSE
356 4 : CALL cp_fm_to_fm(c_z, c_pz)
357 : END IF
358 84 : CALL cp_fm_release(m_tmp)
359 :
360 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo_not_converged, ncol_global=nmo_not_converged, &
361 : context=mo_coeff%matrix_struct%context, &
362 84 : para_env=mo_coeff%matrix_struct%para_env)
363 :
364 84 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm")
365 84 : CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm")
366 84 : CALL cp_fm_struct_release(fm_struct_tmp)
367 :
368 84 : nmat = nmo_not_converged
369 84 : nmat2 = 2*nmo_not_converged
370 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
371 : context=mo_coeff%matrix_struct%context, &
372 84 : para_env=mo_coeff%matrix_struct%para_env)
373 :
374 84 : CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
375 84 : CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
376 84 : CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
377 84 : CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
378 252 : ALLOCATE (evals(nmat2))
379 :
380 84 : CALL cp_fm_struct_release(fm_struct_tmp)
381 :
382 : ! compute CSC
383 84 : CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
384 :
385 : ! compute CHC
386 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp)
387 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1)
388 :
389 : ! compute ZSC
390 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
391 84 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
392 84 : CALL cp_fm_transpose(m_tmp, mt_tmp)
393 84 : CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
394 : ! compute ZHC
395 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
396 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
397 84 : CALL cp_fm_transpose(m_tmp, mt_tmp)
398 84 : CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
399 :
400 84 : CALL cp_fm_release(mt_tmp)
401 :
402 : ! reuse m_sc and m_hc to computr HZ and SZ
403 84 : IF (nmo_converged > 0) THEN
404 12 : CALL parallel_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc)
405 12 : CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc)
406 :
407 12 : CALL cp_fm_release(c_out)
408 12 : CALL cp_fm_release(c_conv)
409 12 : CALL cp_fm_release(s_fm)
410 : ELSE
411 72 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo)
412 72 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo)
413 : END IF
414 :
415 : ! compute ZSZ
416 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
417 84 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
418 : ! compute ZHZ
419 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
420 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
421 :
422 84 : CALL cp_fm_release(m_sc)
423 :
424 : ! solution of the reduced eigenvalues problem
425 84 : CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
426 :
427 : ! extract egenvectors
428 84 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1)
429 84 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc)
430 84 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1)
431 84 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc)
432 :
433 84 : CALL cp_fm_release(m_tmp)
434 :
435 84 : CALL cp_fm_release(c_notconv)
436 84 : CALL cp_fm_release(s_block)
437 84 : CALL cp_fm_release(h_block)
438 84 : CALL cp_fm_release(w_block)
439 84 : CALL cp_fm_release(v_block)
440 :
441 84 : IF (nmo_converged > 0) THEN
442 12 : CALL cp_fm_release(c_z)
443 12 : CALL cp_fm_release(c_pz)
444 12 : DEALLOCATE (c_z, c_pz)
445 12 : jj = 1
446 42 : DO j = 1, nset_not_conv
447 30 : i_first = inotconv_set(j, 1)
448 30 : i_last = inotconv_set(j, 2)
449 30 : n = i_last - i_first + 1
450 30 : CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first)
451 808 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
452 42 : jj = jj + n
453 : END DO
454 12 : DEALLOCATE (iconv_set)
455 12 : DEALLOCATE (inotconv_set)
456 : ELSE
457 72 : CALL cp_fm_to_fm(m_hc, mo_coeff)
458 5328 : eigenvalues(1:nmo) = evals(1:nmo)
459 : END IF
460 84 : DEALLOCATE (evals)
461 :
462 84 : CALL cp_fm_release(m_hc)
463 :
464 84 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
465 :
466 84 : t2 = m_walltime()
467 84 : IF (output_unit > 0) THEN
468 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
469 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
470 : END IF
471 720 : t1 = m_walltime()
472 :
473 : END DO ! iter
474 :
475 40 : DEALLOCATE (iconv)
476 40 : DEALLOCATE (inotconv)
477 40 : DEALLOCATE (ritz_coeff)
478 40 : DEALLOCATE (vnorm)
479 :
480 40 : CALL timestop(handle)
481 120 : END SUBROUTINE generate_extended_space
482 :
483 : ! **************************************************************************************************
484 : !> \brief iterative diagonalization by the block-Davidson approach for one
485 : !> complex K point; complex counterpart of generate_extended_space
486 : !> \param bdav_env Davidson settings of this channel
487 : !> \param mos MO pair of this K point: mos(1) real part, mos(2) imaginary part
488 : !> \param matrix_h complex Kohn-Sham matrix H(k) in full storage
489 : !> \param matrix_s complex overlap matrix S(k) in full storage
490 : !> \param output_unit unit for the Davidson iteration log
491 : !> \param eps_iter convergence threshold for the residual norm of the occupied
492 : !> MOS of this step
493 : !> \param eps_iter_empty threshold for the unoccupied MOS, identical to
494 : !> eps_iter unless the driver relaxed it with EPS_ADAPT
495 : !> \param preconditioner complex K-point preconditioner
496 : !> \note K-point MOs have no sparse dbcsr copy (mo_coeff_b). This routine does
497 : !> not update mo_coeff_b.
498 : ! **************************************************************************************************
499 1906 : SUBROUTINE generate_extended_space_c(bdav_env, mos, matrix_h, matrix_s, output_unit, &
500 : eps_iter, eps_iter_empty, preconditioner)
501 :
502 : TYPE(davidson_type) :: bdav_env
503 : TYPE(mo_set_type), INTENT(INOUT) :: mos
504 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_h, matrix_s
505 : INTEGER, INTENT(IN) :: output_unit
506 : REAL(KIND=dp), INTENT(IN) :: eps_iter, eps_iter_empty
507 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
508 :
509 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space_c'
510 : REAL(KIND=dp), PARAMETER :: eps_exhaust = 1.0E-12_dp, &
511 : occ_tol = 1.0E-3_dp
512 :
513 : COMPLEX(KIND=dp) :: lambda_c
514 1906 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdiag, scaling
515 1906 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: colbuf, submat_buffer
516 : INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, max_iter, n, nao, ncol_z, nmat, &
517 : nmat2, nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
518 1906 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
519 1906 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
520 : LOGICAL :: converged, do_apply_preconditioner, &
521 : exhausted
522 : REAL(KIND=dp) :: eps_iter_col, lambda, max_norm, &
523 : min_norm, t1, t2
524 1906 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm, vnorm_pz
525 1906 : REAL(KIND=dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals, &
526 1906 : occupation
527 : TYPE(cp_cfm_type) :: c_conv, c_mz, c_new, c_notconv, c_out, &
528 : c_pzo, c_pzw, c_z, h_block, h_fm, &
529 : m_hc, m_sc, m_tmp, s_block, v_block, &
530 : w_block
531 : TYPE(cp_cfm_type), POINTER :: cmo
532 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
533 :
534 : ! an occupation below this carries no weight in the density (it sits
535 : ! beyond the thermal smearing tail) and may converge at the relaxed
536 : ! empty threshold
537 :
538 : ! empirical floor for declaring a correction vector numerically zero: the
539 : ! vectors are normalized to one, double-precision roundoff sits near 1e-16,
540 : ! and the threshold errs towards bailing out early (the fallback is an
541 : ! exact diagonalization). Not taken from any reference implementation.
542 :
543 1906 : CALL timeset(routineN, handle)
544 :
545 1906 : NULLIFY (eigenvalues, evals, eig_not_conv, cmo, occupation, fm_struct_tmp)
546 :
547 1906 : do_apply_preconditioner = .FALSE.
548 1906 : IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE.
549 : CALL get_mo_set(mo_set=mos, cmo_coeff=cmo, eigenvalues=eigenvalues, &
550 1906 : occupation_numbers=occupation, nao=nao, nmo=nmo, homo=homo)
551 1906 : IF (do_apply_preconditioner) THEN
552 1778 : max_iter = bdav_env%max_iter
553 : ELSE
554 : max_iter = 1
555 : END IF
556 :
557 1906 : t1 = m_walltime()
558 1906 : IF (output_unit > 0) THEN
559 : WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
560 0 : " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
561 : END IF
562 :
563 5718 : ALLOCATE (iconv(nmo))
564 3812 : ALLOCATE (inotconv(nmo))
565 5718 : ALLOCATE (ritz_coeff(nmo))
566 3812 : ALLOCATE (vnorm(nmo))
567 5718 : ALLOCATE (cdiag(nmo))
568 :
569 1906 : converged = .FALSE.
570 1906 : exhausted = .FALSE.
571 14353 : DO iter = 1, max_iter
572 :
573 : ! compute Ritz values
574 13875 : ritz_coeff = 0.0_dp
575 13875 : CALL cp_cfm_create(m_hc, cmo%matrix_struct, name="hc")
576 13875 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, matrix_h, cmo, z_zero, m_hc)
577 13875 : CALL cp_cfm_create(m_sc, cmo%matrix_struct, name="sc")
578 13875 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, matrix_s, cmo, z_zero, m_sc)
579 :
580 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
581 : context=cmo%matrix_struct%context, &
582 13875 : para_env=cmo%matrix_struct%para_env)
583 13875 : CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
584 13875 : CALL cp_fm_struct_release(fm_struct_tmp)
585 :
586 13875 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmo, m_hc, z_zero, m_tmp)
587 13875 : CALL cp_cfm_get_diag(m_tmp, cdiag)
588 215406 : ritz_coeff(1:nmo) = REAL(cdiag(1:nmo), KIND=dp)
589 13875 : CALL cp_cfm_release(m_tmp)
590 :
591 : ! Check for converged eigenvectors: residual C_z = S*C*diag(ritz) - H*C
592 13875 : CALL cp_cfm_create(c_z, cmo%matrix_struct, name="z")
593 13875 : CALL cp_cfm_to_cfm(m_sc, c_z)
594 13875 : IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
595 27750 : ALLOCATE (scaling(nmo))
596 215406 : scaling(:) = CMPLX(ritz_coeff, 0.0_dp, KIND=dp)
597 13875 : CALL cp_cfm_column_scale(c_z, scaling)
598 13875 : CALL cp_cfm_scale_and_add(-z_one, c_z, z_one, m_hc)
599 13875 : CALL cp_cfm_vectorsnorm(c_z, vnorm)
600 :
601 13875 : nmo_converged = 0
602 13875 : nmo_not_converged = 0
603 13875 : max_norm = 0.0_dp
604 13875 : min_norm = 1.e10_dp
605 215406 : DO imo = 1, nmo
606 201531 : max_norm = MAX(max_norm, vnorm(imo))
607 215406 : min_norm = MIN(min_norm, vnorm(imo))
608 : END DO
609 13875 : iconv = 0
610 13875 : inotconv = 0
611 215406 : DO imo = 1, nmo
612 : ! the relaxed threshold applies to the UNOCCUPIED manifold only:
613 : ! under smearing the MOS above the HOMO can still carry weight and
614 : ! leak their residual into the density, so the occupation decides,
615 : ! not the HOMO index
616 201531 : IF (occupation(imo) < occ_tol) THEN
617 126420 : eps_iter_col = eps_iter_empty
618 : ELSE
619 75111 : eps_iter_col = eps_iter
620 : END IF
621 215406 : IF (vnorm(imo) <= eps_iter_col) THEN
622 4690 : nmo_converged = nmo_converged + 1
623 4690 : iconv(nmo_converged) = imo
624 : ELSE
625 196841 : nmo_not_converged = nmo_not_converged + 1
626 196841 : inotconv(nmo_not_converged) = imo
627 : END IF
628 : END DO
629 :
630 : ! the iter > 1 gate below blocks the exit when EVERY column passes
631 : ! the entering check at iter == 1. With no unconverged column there
632 : ! is no correction to rotate with. The reduced problem further down
633 : ! would be zero-dimensional, and this head's nmo-wide residual
634 : ! matrices would be released for nothing by the packing block below.
635 : ! Reclassify all columns as unconverged so the call still performs
636 : ! the required one Rayleigh-Ritz rotation under the current
637 : ! operator: the zero-converged branch reuses the head matrices
638 : ! directly. At iter > 1 the same situation exits legitimately above,
639 : ! hence the iter guard
640 13875 : IF (iter == 1 .AND. nmo_not_converged == 0) THEN
641 244 : nmo_not_converged = nmo
642 244 : nmo_converged = 0
643 244 : iconv = 0
644 492 : DO imo = 1, nmo
645 492 : inotconv(imo) = imo
646 : END DO
647 : END IF
648 :
649 13875 : IF (nmo_converged > 0) THEN
650 5292 : ALLOCATE (iconv_set(nmo_converged, 2))
651 4040 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
652 1764 : i_last = iconv(1)
653 1764 : nset = 0
654 6206 : DO j = 1, nmo_converged
655 4442 : imo = iconv(j)
656 :
657 6206 : IF (imo == i_last + 1) THEN
658 2008 : i_last = imo
659 2008 : iconv_set(nset, 2) = imo
660 : ELSE
661 2434 : i_last = imo
662 2434 : nset = nset + 1
663 2434 : iconv_set(nset, 1) = imo
664 2434 : iconv_set(nset, 2) = imo
665 : END IF
666 : END DO
667 1764 : nset_conv = nset
668 :
669 1764 : i_last = inotconv(1)
670 1764 : nset = 0
671 7858 : DO j = 1, nmo_not_converged
672 6094 : imo = inotconv(j)
673 :
674 7858 : IF (imo == i_last + 1) THEN
675 4900 : i_last = imo
676 4900 : inotconv_set(nset, 2) = imo
677 : ELSE
678 1194 : i_last = imo
679 1194 : nset = nset + 1
680 1194 : inotconv_set(nset, 1) = imo
681 1194 : inotconv_set(nset, 2) = imo
682 : END IF
683 : END DO
684 1764 : nset_not_conv = nset
685 1764 : CALL cp_cfm_release(m_sc)
686 1764 : CALL cp_cfm_release(m_hc)
687 1764 : CALL cp_cfm_release(c_z)
688 : END IF
689 :
690 : ! the convergence check runs at the head of a cycle and probes the
691 : ! ENTERING MOS. Requiring iter > 1 enforces one Rayleigh-Ritz rotation
692 : ! under the current operator per call. Exiting on the entry check
693 : ! alone would return the MOS unchanged, the density unchanged, and an
694 : ! outer SCF convergence measure of exactly zero. The SCF then stops
695 : ! at whatever the previous step left, which masquerades as
696 : ! convergence whenever the threshold is looser than the entering
697 : ! residual
698 13875 : IF (iter > 1 .AND. REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
699 1280 : converged = .TRUE.
700 1280 : DEALLOCATE (iconv_set)
701 1280 : DEALLOCATE (inotconv_set)
702 1280 : t2 = m_walltime()
703 1280 : IF (output_unit > 0) THEN
704 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
705 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
706 :
707 0 : WRITE (output_unit, *) " Reached convergence in ", iter, &
708 0 : " Davidson iterations"
709 : END IF
710 :
711 : EXIT
712 : END IF
713 :
714 12595 : ncol_z = nmo_not_converged
715 :
716 12595 : IF (nmo_converged > 0) THEN
717 : ! dense copy of H, shifted below by the projector onto the converged space
718 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
719 : context=cmo%matrix_struct%context, &
720 484 : para_env=cmo%matrix_struct%para_env)
721 484 : CALL cp_cfm_create(h_fm, fm_struct_tmp, name="hf")
722 484 : CALL cp_cfm_create(c_out, fm_struct_tmp, name="cout")
723 484 : CALL cp_fm_struct_release(fm_struct_tmp)
724 484 : CALL cp_cfm_to_cfm(matrix_h, h_fm)
725 :
726 : ! projector P = (S*C_conv)*(S*C_conv)^H
727 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
728 : context=cmo%matrix_struct%context, &
729 484 : para_env=cmo%matrix_struct%para_env)
730 484 : CALL cp_cfm_create(c_conv, fm_struct_tmp, name="cconv")
731 484 : CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="scconv")
732 484 : CALL cp_fm_struct_release(fm_struct_tmp)
733 484 : jj = 1
734 1618 : DO j = 1, nset_conv
735 1134 : i_first = iconv_set(j, 1)
736 1134 : i_last = iconv_set(j, 2)
737 1134 : n = i_last - i_first + 1
738 4536 : ALLOCATE (submat_buffer(nao, n))
739 : CALL cp_cfm_get_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first, &
740 1134 : n_rows=nao, n_cols=n)
741 1134 : CALL cp_cfm_set_submatrix(c_conv, submat_buffer, start_row=1, start_col=jj)
742 1134 : DEALLOCATE (submat_buffer)
743 1618 : jj = jj + n
744 : END DO
745 : CALL cp_cfm_gemm('N', 'N', nao, nmo_converged, nao, z_one, matrix_s, c_conv, &
746 484 : z_zero, m_tmp)
747 484 : CALL cp_cfm_gemm('N', 'C', nao, nao, nmo_converged, z_one, m_tmp, m_tmp, z_zero, c_out)
748 484 : CALL cp_cfm_release(m_tmp)
749 484 : CALL cp_cfm_release(c_conv)
750 :
751 : ! project c_out out of H
752 484 : lambda = 100.0_dp*ABS(eigenvalues(homo))
753 484 : lambda_c = CMPLX(lambda, 0.0_dp, KIND=dp)
754 484 : CALL cp_cfm_scale_and_add(z_one, h_fm, lambda_c, c_out)
755 : END IF
756 :
757 : ! gather the not converged MOs
758 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=ncol_z, &
759 : context=cmo%matrix_struct%context, &
760 12595 : para_env=cmo%matrix_struct%para_env)
761 12595 : CALL cp_cfm_create(c_notconv, fm_struct_tmp, name="c_notconv")
762 12595 : CALL cp_fm_struct_release(fm_struct_tmp)
763 :
764 12595 : IF (nmo_converged > 0) THEN
765 : ! m_hc/m_sc/c_z are recreated with the ncol_z column count
766 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=ncol_z, &
767 : context=cmo%matrix_struct%context, &
768 484 : para_env=cmo%matrix_struct%para_env)
769 484 : CALL cp_cfm_create(m_hc, fm_struct_tmp, name="m_hc")
770 484 : CALL cp_cfm_create(m_sc, fm_struct_tmp, name="m_sc")
771 484 : CALL cp_cfm_create(c_z, fm_struct_tmp, name="c_z")
772 484 : CALL cp_fm_struct_release(fm_struct_tmp)
773 :
774 1452 : ALLOCATE (eig_not_conv(ncol_z))
775 484 : jj = 1
776 1634 : DO j = 1, nset_not_conv
777 1150 : i_first = inotconv_set(j, 1)
778 1150 : i_last = inotconv_set(j, 2)
779 1150 : n = i_last - i_first + 1
780 4600 : ALLOCATE (submat_buffer(nao, n))
781 : CALL cp_cfm_get_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first, &
782 1150 : n_rows=nao, n_cols=n)
783 1150 : CALL cp_cfm_set_submatrix(c_notconv, submat_buffer, start_row=1, start_col=jj)
784 1150 : DEALLOCATE (submat_buffer)
785 7172 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
786 1634 : jj = jj + n
787 : END DO
788 : ! extend the subspace using only the not converged vectors
789 484 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, h_fm, c_notconv, z_zero, m_hc)
790 484 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_s, c_notconv, z_zero, m_sc)
791 484 : CALL cp_cfm_to_cfm(m_sc, c_z)
792 484 : IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
793 1452 : ALLOCATE (scaling(ncol_z))
794 6506 : scaling(:) = CMPLX(eig_not_conv, 0.0_dp, KIND=dp)
795 484 : CALL cp_cfm_column_scale(c_z, scaling)
796 484 : CALL cp_cfm_scale_and_add(-z_one, c_z, z_one, m_hc)
797 968 : DEALLOCATE (eig_not_conv)
798 : ! h_fm is the shifted operator H+lambda*P, still needed for the HZ block below
799 : ELSE
800 : ! nothing frozen: c_notconv is the full set and the m_hc/m_sc/c_z matrices
801 : ! of the Ritz step above (nmo-wide) are reused directly
802 12111 : CALL cp_cfm_to_cfm(cmo, c_notconv)
803 : END IF
804 :
805 : ! preconditioner
806 12595 : CALL cp_cfm_create(c_mz, c_z%matrix_struct, name="pz")
807 12595 : IF (do_apply_preconditioner) THEN
808 12467 : IF (preconditioner%in_use /= 0) THEN
809 12467 : IF (nmo_converged == 0) THEN
810 : ! every column is unconverged, so the residual set already is
811 : ! the nmo-wide input the applier contract requires: apply
812 : ! directly, no scattering
813 11983 : CALL apply_preconditioner(preconditioner, c_z, c_mz)
814 : ELSE
815 : ! the complex applier works column-wise in the MO basis of
816 : ! the preconditioner construction, so the unconverged
817 : ! residuals are scattered into an nmo-wide zero buffer and
818 : ! gathered back. The unconverged columns move per contiguous
819 : ! segment of inotconv_set. The submatrix transfers end in a
820 : ! collective reduction over the matrix group, and a
821 : ! column-at-a-time loop would launch one per unconverged MO
822 484 : CALL cp_cfm_create(c_pzw, cmo%matrix_struct, name="pzw")
823 484 : CALL cp_cfm_create(c_pzo, cmo%matrix_struct, name="pzo")
824 484 : CALL cp_cfm_set_all(c_pzw, z_zero, z_zero)
825 1936 : ALLOCATE (colbuf(nao, ncol_z))
826 484 : jj = 1
827 1634 : DO j = 1, nset_not_conv
828 1150 : i_first = inotconv_set(j, 1)
829 1150 : i_last = inotconv_set(j, 2)
830 1150 : n = i_last - i_first + 1
831 1150 : CALL cp_cfm_get_submatrix(c_z, colbuf(:, 1:n), n_rows=nao, n_cols=n, start_col=jj)
832 : CALL cp_cfm_set_submatrix(c_pzw, colbuf(:, 1:n), n_rows=nao, n_cols=n, &
833 1150 : start_col=i_first)
834 1634 : jj = jj + n
835 : END DO
836 484 : CALL apply_preconditioner(preconditioner, c_pzw, c_pzo)
837 484 : jj = 1
838 1634 : DO j = 1, nset_not_conv
839 1150 : i_first = inotconv_set(j, 1)
840 1150 : i_last = inotconv_set(j, 2)
841 1150 : n = i_last - i_first + 1
842 : CALL cp_cfm_get_submatrix(c_pzo, colbuf(:, 1:n), n_rows=nao, n_cols=n, &
843 1150 : start_col=i_first)
844 1150 : CALL cp_cfm_set_submatrix(c_mz, colbuf(:, 1:n), n_rows=nao, n_cols=n, start_col=jj)
845 1634 : jj = jj + n
846 : END DO
847 484 : DEALLOCATE (colbuf)
848 484 : CALL cp_cfm_release(c_pzw)
849 968 : CALL cp_cfm_release(c_pzo)
850 : END IF
851 : ELSE
852 0 : CALL cp_cfm_to_cfm(c_z, c_mz)
853 : END IF
854 : ELSE
855 128 : CALL cp_cfm_to_cfm(c_z, c_mz)
856 : END IF
857 :
858 : ! normalize the correction vectors, then remove their components inside
859 : ! the span of the current MOS, z <- z - C*(C^H*S*z): the blocked Davidson
860 : ! subspace expansion of Kresse and Furthmueller (1996). This keeps the
861 : ! reduced S block well conditioned when a MOS is already (nearly) an
862 : ! eigenvector and its preconditioned residual is (nearly) parallel to it.
863 : ! Exhaustion: either the preconditioned correction itself, or its
864 : ! component outside the MOS span, can vanish by symmetry (the density
865 : ! keeps the little-group symmetry of the kpoint), which would make the
866 : ! reduced overlap block singular. Stop the iteration and finish the
867 : ! channel by the direct diagonalization below. The scratch of the
868 : ! interrupted iteration is released at the single cleanup point after
869 : ! the loop.
870 37785 : ALLOCATE (vnorm_pz(ncol_z))
871 12595 : CALL cp_cfm_vectorsnorm(c_mz, vnorm_pz)
872 209612 : IF (ANY(vnorm_pz < eps_exhaust)) THEN
873 : exhausted = .TRUE.
874 : ELSE
875 12595 : IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
876 37785 : ALLOCATE (scaling(ncol_z))
877 209612 : scaling(:) = CMPLX(1.0_dp/vnorm_pz, 0.0_dp, KIND=dp)
878 12595 : CALL cp_cfm_column_scale(c_mz, scaling)
879 :
880 : ! projection coefficients C^H*S*z. c_z is reused as scratch for S*z
881 12595 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_s, c_mz, z_zero, c_z)
882 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=ncol_z, &
883 : context=cmo%matrix_struct%context, &
884 12595 : para_env=cmo%matrix_struct%para_env)
885 12595 : CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="proj_coef")
886 12595 : CALL cp_fm_struct_release(fm_struct_tmp)
887 12595 : CALL cp_cfm_gemm('C', 'N', nmo, ncol_z, nao, z_one, cmo, c_z, z_zero, m_tmp)
888 12595 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nmo, -z_one, cmo, m_tmp, z_one, c_mz)
889 12595 : CALL cp_cfm_release(m_tmp)
890 :
891 12595 : CALL cp_cfm_vectorsnorm(c_mz, vnorm_pz)
892 221911 : IF (ANY(vnorm_pz < eps_exhaust)) THEN
893 : exhausted = .TRUE.
894 : ELSE
895 12447 : IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
896 24894 : ALLOCATE (scaling(ncol_z))
897 209168 : scaling(:) = CMPLX(1.0_dp/vnorm_pz, 0.0_dp, KIND=dp)
898 12447 : CALL cp_cfm_column_scale(c_mz, scaling)
899 : END IF
900 : END IF
901 12595 : DEALLOCATE (vnorm_pz)
902 12595 : IF (exhausted) EXIT
903 :
904 12447 : nmat = ncol_z
905 12447 : nmat2 = 2*ncol_z
906 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
907 : context=cmo%matrix_struct%context, &
908 12447 : para_env=cmo%matrix_struct%para_env)
909 :
910 12447 : CALL cp_cfm_create(s_block, fm_struct_tmp, name="sb")
911 12447 : CALL cp_cfm_create(h_block, fm_struct_tmp, name="hb")
912 12447 : CALL cp_cfm_create(v_block, fm_struct_tmp, name="vb")
913 12447 : CALL cp_cfm_create(w_block, fm_struct_tmp, name="wb")
914 37341 : ALLOCATE (evals(nmat2))
915 12447 : CALL cp_fm_struct_release(fm_struct_tmp)
916 :
917 : ! compute CHC and CSC first, with m_hc/m_sc still holding the
918 : ! H*C_nc and S*C_nc products of the residual step
919 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_notconv, m_hc, z_zero, h_block, &
920 12447 : c_first_row=1, c_first_col=1)
921 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_notconv, m_sc, z_zero, s_block, &
922 12447 : c_first_row=1, c_first_col=1)
923 :
924 : ! compute ZSC and ZHC with m_sc/m_sc still holding S*C_nc and H*C_nc
925 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_mz, m_sc, z_zero, s_block, &
926 12447 : c_first_row=1 + nmat, c_first_col=1)
927 : ! compute ZHC
928 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_mz, m_hc, z_zero, h_block, &
929 12447 : c_first_row=1 + nmat, c_first_col=1)
930 :
931 : ! then reuse m_sc and m_hc to compute SZ and HZ (mirrors the Gamma
932 : ! version); with frozen converged vectors the operator is the shifted
933 : ! h_fm = H + lambda*P
934 12447 : IF (nmo_converged > 0) THEN
935 484 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, h_fm, c_mz, z_zero, m_hc)
936 484 : CALL cp_cfm_release(c_out)
937 484 : CALL cp_cfm_release(h_fm)
938 : ELSE
939 11963 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_h, c_mz, z_zero, m_hc)
940 : END IF
941 12447 : CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_s, c_mz, z_zero, m_sc)
942 :
943 : ! the opposite off-diagonal blocks are the conjugate transposes
944 : ! (ZSC)^H = C_nc^H*(S*Z) and (ZHC)^H = C_nc^H*(H*Z): with S and H
945 : ! Hermitian these are the same products read in the other direction,
946 : ! so the blocks complete as gemms against the just-computed S*Z and
947 : ! H*Z, with no gather/conjugate/scatter round trip per block (two
948 : ! collectives and a serial conjugation of nmat^2 elements each)
949 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_notconv, m_sc, z_zero, s_block, &
950 12447 : c_first_row=1, c_first_col=1 + nmat)
951 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_notconv, m_hc, z_zero, h_block, &
952 12447 : c_first_row=1, c_first_col=1 + nmat)
953 :
954 : ! compute ZSZ
955 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_mz, m_sc, z_zero, s_block, &
956 12447 : c_first_row=1 + nmat, c_first_col=1 + nmat)
957 : ! compute ZHZ
958 : CALL cp_cfm_gemm('C', 'N', nmat, nmat, nao, z_one, c_mz, m_hc, z_zero, h_block, &
959 12447 : c_first_row=1 + nmat, c_first_col=1 + nmat)
960 :
961 12447 : CALL cp_cfm_release(m_sc)
962 :
963 : ! solution of the reduced generalized eigenproblem
964 12447 : CALL cp_cfm_geeig(h_block, s_block, v_block, evals, w_block)
965 :
966 : ! new MOS: C_nc*V_11 + Z*V_21
967 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmat, &
968 : context=cmo%matrix_struct%context, &
969 12447 : para_env=cmo%matrix_struct%para_env)
970 12447 : CALL cp_cfm_create(c_new, fm_struct_tmp, name="c_new")
971 12447 : CALL cp_fm_struct_release(fm_struct_tmp)
972 : CALL cp_cfm_gemm('N', 'N', nao, nmat, nmat, z_one, c_notconv, v_block, z_zero, c_new, &
973 12447 : b_first_row=1, b_first_col=1)
974 : CALL cp_cfm_gemm('N', 'N', nao, nmat, nmat, z_one, c_mz, v_block, z_one, c_new, &
975 12447 : b_first_row=1 + nmat, b_first_col=1)
976 :
977 12447 : CALL cp_cfm_release(m_hc)
978 12447 : CALL cp_cfm_release(c_mz)
979 12447 : CALL cp_cfm_release(c_z)
980 12447 : CALL cp_cfm_release(c_notconv)
981 12447 : CALL cp_cfm_release(s_block)
982 12447 : CALL cp_cfm_release(h_block)
983 12447 : CALL cp_cfm_release(w_block)
984 12447 : CALL cp_cfm_release(v_block)
985 :
986 12447 : IF (nmo_converged > 0) THEN
987 484 : jj = 1
988 1634 : DO j = 1, nset_not_conv
989 1150 : i_first = inotconv_set(j, 1)
990 1150 : i_last = inotconv_set(j, 2)
991 1150 : n = i_last - i_first + 1
992 4600 : ALLOCATE (submat_buffer(nao, n))
993 : CALL cp_cfm_get_submatrix(c_new, submat_buffer, start_row=1, start_col=jj, &
994 1150 : n_rows=nao, n_cols=n)
995 1150 : CALL cp_cfm_set_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first)
996 1150 : DEALLOCATE (submat_buffer)
997 14344 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
998 1634 : jj = jj + n
999 : END DO
1000 484 : DEALLOCATE (iconv_set)
1001 484 : DEALLOCATE (inotconv_set)
1002 : ELSE
1003 11963 : CALL cp_cfm_to_cfm(c_new, cmo)
1004 405324 : eigenvalues(1:nmo) = evals(1:nmo)
1005 : END IF
1006 12447 : DEALLOCATE (evals)
1007 12447 : CALL cp_cfm_release(c_new)
1008 :
1009 12447 : t2 = m_walltime()
1010 12447 : IF (output_unit > 0) THEN
1011 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1012 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
1013 : END IF
1014 118657 : t1 = m_walltime()
1015 :
1016 : END DO ! iter
1017 :
1018 1906 : IF (exhausted) THEN
1019 : ! single cleanup point for the interrupted iteration: the shifted
1020 : ! operator h_fm ends its life either here or at the HZ recompute
1021 : ! inside the loop, its two end-of-life sites
1022 148 : IF (nmo_converged > 0) THEN
1023 0 : CALL cp_cfm_release(h_fm)
1024 0 : CALL cp_cfm_release(c_out)
1025 : END IF
1026 148 : CALL cp_cfm_release(m_hc)
1027 148 : CALL cp_cfm_release(m_sc)
1028 148 : CALL cp_cfm_release(c_z)
1029 148 : CALL cp_cfm_release(c_mz)
1030 148 : CALL cp_cfm_release(c_notconv)
1031 148 : IF (ALLOCATED(iconv_set)) DEALLOCATE (iconv_set)
1032 148 : IF (ALLOCATED(inotconv_set)) DEALLOCATE (inotconv_set)
1033 :
1034 : ! the subspace expansion was exhausted (symmetry-invariant MOS span);
1035 : ! finish this channel by a direct diagonalization of the intact
1036 : ! operator, so that at least the rotation inside the MOS span is exact
1037 148 : CALL cp_cfm_create(h_fm, matrix_h%matrix_struct, name="fb_h")
1038 148 : CALL cp_cfm_create(c_out, matrix_h%matrix_struct, name="fb_s")
1039 148 : CALL cp_cfm_create(c_new, matrix_h%matrix_struct, name="fb_work")
1040 148 : CALL cp_cfm_to_cfm(matrix_h, h_fm)
1041 148 : CALL cp_cfm_to_cfm(matrix_s, c_out)
1042 148 : CALL cp_cfm_geeig(h_fm, c_out, cmo, eigenvalues, c_new)
1043 148 : CALL cp_cfm_release(h_fm)
1044 148 : CALL cp_cfm_release(c_out)
1045 148 : CALL cp_cfm_release(c_new)
1046 : END IF
1047 :
1048 1906 : DEALLOCATE (iconv)
1049 1906 : DEALLOCATE (inotconv)
1050 1906 : DEALLOCATE (ritz_coeff)
1051 1906 : DEALLOCATE (vnorm)
1052 1906 : DEALLOCATE (cdiag)
1053 1906 : IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
1054 :
1055 1906 : CALL timestop(handle)
1056 5718 : END SUBROUTINE generate_extended_space_c
1057 :
1058 : ! **************************************************************************************************
1059 : !> \brief ...
1060 : !> \param bdav_env ...
1061 : !> \param mo_set ...
1062 : !> \param matrix_h ...
1063 : !> \param matrix_s ...
1064 : !> \param output_unit ...
1065 : !> \param preconditioner ...
1066 : ! **************************************************************************************************
1067 64 : SUBROUTINE generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
1068 : preconditioner)
1069 :
1070 : TYPE(davidson_type) :: bdav_env
1071 : TYPE(mo_set_type), INTENT(IN) :: mo_set
1072 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
1073 : INTEGER, INTENT(IN) :: output_unit
1074 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
1075 :
1076 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space_sparse'
1077 :
1078 : INTEGER :: col_offset, handle, homo, i_first, i_last, imo, iteration, j, jj, k, max_iter, n, &
1079 : nao, nmat, nmat2, nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
1080 64 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
1081 64 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
1082 : LOGICAL :: converged, do_apply_preconditioner
1083 : REAL(dp) :: lambda, max_norm, min_norm, t1, t2
1084 64 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, ritz_coeff, vnorm
1085 64 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues
1086 64 : REAL(dp), DIMENSION(:, :), POINTER :: block
1087 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1088 : TYPE(cp_fm_type) :: h_block, matrix_mm_fm, matrix_mmt_fm, &
1089 : matrix_nm_fm, matrix_z_fm, mo_conv_fm, &
1090 : s_block, v_block, w_block
1091 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_notconv_fm
1092 : TYPE(dbcsr_iterator_type) :: iter
1093 : TYPE(dbcsr_type), POINTER :: c_out, matrix_hc, matrix_mm, matrix_pz, &
1094 : matrix_sc, matrix_z, mo_coeff_b, &
1095 : mo_conv, mo_notconv, smo_conv
1096 : TYPE(mp_comm_type) :: group
1097 :
1098 64 : CALL timeset(routineN, handle)
1099 :
1100 64 : do_apply_preconditioner = .FALSE.
1101 64 : IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE.
1102 :
1103 64 : NULLIFY (mo_coeff, mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm)
1104 64 : NULLIFY (mo_notconv_fm, mo_conv, mo_notconv, smo_conv, c_out)
1105 64 : NULLIFY (fm_struct_tmp)
1106 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
1107 64 : eigenvalues=eigenvalues, homo=homo, nao=nao, nmo=nmo)
1108 64 : IF (do_apply_preconditioner) THEN
1109 56 : max_iter = bdav_env%max_iter
1110 : ELSE
1111 : max_iter = 1
1112 : END IF
1113 :
1114 64 : t1 = m_walltime()
1115 64 : IF (output_unit > 0) THEN
1116 : WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
1117 0 : " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
1118 : END IF
1119 :
1120 : ! Allocate array for Ritz values
1121 192 : ALLOCATE (ritz_coeff(nmo))
1122 192 : ALLOCATE (iconv(nmo))
1123 128 : ALLOCATE (inotconv(nmo))
1124 128 : ALLOCATE (vnorm(nmo))
1125 :
1126 64 : converged = .FALSE.
1127 128 : DO iteration = 1, max_iter
1128 88 : NULLIFY (c_out, mo_conv, mo_notconv_fm, mo_notconv)
1129 : ! Prepare HC and SC, using mo_coeff_b (sparse), these are still sparse
1130 88 : CALL dbcsr_init_p(matrix_hc)
1131 : CALL dbcsr_create(matrix_hc, template=mo_coeff_b, &
1132 : name="matrix_hc", &
1133 88 : matrix_type=dbcsr_type_no_symmetry)
1134 88 : CALL dbcsr_init_p(matrix_sc)
1135 : CALL dbcsr_create(matrix_sc, template=mo_coeff_b, &
1136 : name="matrix_sc", &
1137 88 : matrix_type=dbcsr_type_no_symmetry)
1138 :
1139 88 : CALL dbcsr_get_info(mo_coeff_b, nfullrows_total=n, nfullcols_total=k, group=group)
1140 88 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h, mo_coeff_b, 0.0_dp, matrix_hc, last_column=k)
1141 88 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_coeff_b, 0.0_dp, matrix_sc, last_column=k)
1142 :
1143 : ! compute Ritz values
1144 88 : ritz_coeff = 0.0_dp
1145 : ! Allocate Sparse matrices: nmoxnmo
1146 : ! matrix_mm
1147 :
1148 88 : CALL dbcsr_init_p(matrix_mm)
1149 : CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo, n=nmo, &
1150 88 : sym=dbcsr_type_no_symmetry)
1151 :
1152 88 : CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_coeff_b, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
1153 88 : CALL dbcsr_get_diag(matrix_mm, ritz_coeff)
1154 88 : CALL mo_coeff%matrix_struct%para_env%sum(ritz_coeff)
1155 :
1156 : ! extended subspace P Z = P [H - theta S]C this ia another matrix of type and size as mo_coeff_b
1157 88 : CALL dbcsr_init_p(matrix_z)
1158 : CALL dbcsr_create(matrix_z, template=mo_coeff_b, &
1159 : name="matrix_z", &
1160 88 : matrix_type=dbcsr_type_no_symmetry)
1161 88 : CALL dbcsr_copy(matrix_z, matrix_sc)
1162 88 : CALL dbcsr_scale_by_vector(matrix_z, ritz_coeff, side='right')
1163 88 : CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
1164 :
1165 : ! Compute the column norms of matrix_z.
1166 88 : vnorm = 0.0_dp
1167 88 : CALL dbcsr_iterator_start(iter, matrix_z)
1168 792 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1169 704 : CALL dbcsr_iterator_next_block(iter, block=block, col_offset=col_offset)
1170 13464 : DO j = 1, SIZE(block, 2)
1171 178112 : vnorm(col_offset + j - 1) = vnorm(col_offset + j - 1) + SUM(block(:, j)**2)
1172 : END DO
1173 : END DO
1174 88 : CALL dbcsr_iterator_stop(iter)
1175 88 : CALL group%sum(vnorm)
1176 3256 : vnorm = SQRT(vnorm)
1177 :
1178 : ! Check for converged eigenvectors
1179 88 : nmo_converged = 0
1180 88 : nmo_not_converged = 0
1181 88 : max_norm = 0.0_dp
1182 88 : min_norm = 1.e10_dp
1183 3256 : DO imo = 1, nmo
1184 3168 : max_norm = MAX(max_norm, vnorm(imo))
1185 3256 : min_norm = MIN(min_norm, vnorm(imo))
1186 : END DO
1187 88 : iconv = 0
1188 88 : inotconv = 0
1189 :
1190 3256 : DO imo = 1, nmo
1191 3256 : IF (vnorm(imo) <= bdav_env%eps_iter) THEN
1192 836 : nmo_converged = nmo_converged + 1
1193 836 : iconv(nmo_converged) = imo
1194 : ELSE
1195 2332 : nmo_not_converged = nmo_not_converged + 1
1196 2332 : inotconv(nmo_not_converged) = imo
1197 : END IF
1198 : END DO
1199 :
1200 88 : IF (nmo_converged > 0) THEN
1201 90 : ALLOCATE (iconv_set(nmo_converged, 2))
1202 88 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
1203 30 : i_last = iconv(1)
1204 30 : nset = 0
1205 866 : DO j = 1, nmo_converged
1206 836 : imo = iconv(j)
1207 :
1208 866 : IF (imo == i_last + 1) THEN
1209 772 : i_last = imo
1210 772 : iconv_set(nset, 2) = imo
1211 : ELSE
1212 64 : i_last = imo
1213 64 : nset = nset + 1
1214 64 : iconv_set(nset, 1) = imo
1215 64 : iconv_set(nset, 2) = imo
1216 : END IF
1217 : END DO
1218 30 : nset_conv = nset
1219 :
1220 30 : i_last = inotconv(1)
1221 30 : nset = 0
1222 274 : DO j = 1, nmo_not_converged
1223 244 : imo = inotconv(j)
1224 :
1225 274 : IF (imo == i_last + 1) THEN
1226 184 : i_last = imo
1227 184 : inotconv_set(nset, 2) = imo
1228 : ELSE
1229 60 : i_last = imo
1230 60 : nset = nset + 1
1231 60 : inotconv_set(nset, 1) = imo
1232 60 : inotconv_set(nset, 2) = imo
1233 : END IF
1234 : END DO
1235 30 : nset_not_conv = nset
1236 :
1237 30 : CALL dbcsr_release_p(matrix_hc)
1238 30 : CALL dbcsr_release_p(matrix_sc)
1239 30 : CALL dbcsr_release_p(matrix_z)
1240 30 : CALL dbcsr_release_p(matrix_mm)
1241 : END IF
1242 :
1243 88 : IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
1244 24 : DEALLOCATE (iconv_set)
1245 :
1246 24 : DEALLOCATE (inotconv_set)
1247 :
1248 24 : converged = .TRUE.
1249 24 : t2 = m_walltime()
1250 24 : IF (output_unit > 0) THEN
1251 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1252 0 : iteration, nmo_converged, max_norm, min_norm, t2 - t1
1253 :
1254 0 : WRITE (output_unit, *) " Reached convergence in ", iteration, &
1255 0 : " Davidson iterations"
1256 : END IF
1257 :
1258 : EXIT
1259 : END IF
1260 :
1261 64 : IF (nmo_converged > 0) THEN
1262 :
1263 : !allocate mo_conv_fm
1264 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
1265 : context=mo_coeff%matrix_struct%context, &
1266 6 : para_env=mo_coeff%matrix_struct%para_env)
1267 6 : CALL cp_fm_create(mo_conv_fm, fm_struct_tmp, name="mo_conv_fm")
1268 :
1269 6 : CALL cp_fm_struct_release(fm_struct_tmp)
1270 :
1271 : ! extract mo_conv from mo_coeff full matrix
1272 6 : jj = 1
1273 22 : DO j = 1, nset_conv
1274 16 : i_first = iconv_set(j, 1)
1275 16 : i_last = iconv_set(j, 2)
1276 16 : n = i_last - i_first + 1
1277 16 : CALL cp_fm_to_fm_submat(mo_coeff, mo_conv_fm, nao, n, 1, i_first, 1, jj)
1278 22 : jj = jj + n
1279 : END DO
1280 :
1281 : ! allocate c_out sparse matrix, to project out the converged MOS
1282 6 : CALL dbcsr_init_p(c_out)
1283 : CALL dbcsr_create(c_out, template=matrix_s, &
1284 : name="c_out", &
1285 6 : matrix_type=dbcsr_type_symmetric)
1286 :
1287 : ! allocate mo_conv sparse
1288 6 : CALL dbcsr_init_p(mo_conv)
1289 : CALL cp_dbcsr_m_by_n_from_row_template(mo_conv, template=matrix_s, n=nmo_converged, &
1290 6 : sym=dbcsr_type_no_symmetry)
1291 :
1292 6 : CALL dbcsr_init_p(smo_conv)
1293 : CALL cp_dbcsr_m_by_n_from_row_template(smo_conv, template=matrix_s, n=nmo_converged, &
1294 6 : sym=dbcsr_type_no_symmetry)
1295 :
1296 6 : CALL copy_fm_to_dbcsr(mo_conv_fm, mo_conv) !fm->dbcsr
1297 :
1298 6 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_conv, 0.0_dp, smo_conv, last_column=nmo_converged)
1299 6 : CALL dbcsr_multiply('n', 't', 1.0_dp, smo_conv, smo_conv, 0.0_dp, c_out, last_column=nao)
1300 : ! project c_out out of H
1301 6 : lambda = 100.0_dp*ABS(eigenvalues(homo))
1302 6 : CALL dbcsr_add(c_out, matrix_h, lambda, 1.0_dp)
1303 :
1304 6 : CALL dbcsr_release_p(mo_conv)
1305 6 : CALL dbcsr_release_p(smo_conv)
1306 6 : CALL cp_fm_release(mo_conv_fm)
1307 :
1308 : !allocate c_notconv_fm
1309 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
1310 : context=mo_coeff%matrix_struct%context, &
1311 6 : para_env=mo_coeff%matrix_struct%para_env)
1312 6 : ALLOCATE (mo_notconv_fm)
1313 6 : CALL cp_fm_create(mo_notconv_fm, fm_struct_tmp, name="mo_notconv_fm")
1314 6 : CALL cp_fm_struct_release(fm_struct_tmp)
1315 :
1316 : ! extract mo_notconv from mo_coeff full matrix
1317 18 : ALLOCATE (eig_not_conv(nmo_not_converged))
1318 :
1319 6 : jj = 1
1320 24 : DO j = 1, nset_not_conv
1321 18 : i_first = inotconv_set(j, 1)
1322 18 : i_last = inotconv_set(j, 2)
1323 18 : n = i_last - i_first + 1
1324 18 : CALL cp_fm_to_fm_submat(mo_coeff, mo_notconv_fm, nao, n, 1, i_first, 1, jj)
1325 186 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
1326 24 : jj = jj + n
1327 : END DO
1328 :
1329 : ! allocate mo_conv sparse
1330 6 : CALL dbcsr_init_p(mo_notconv)
1331 : CALL cp_dbcsr_m_by_n_from_row_template(mo_notconv, template=matrix_s, n=nmo_not_converged, &
1332 6 : sym=dbcsr_type_no_symmetry)
1333 :
1334 6 : CALL dbcsr_init_p(matrix_hc)
1335 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_hc, template=matrix_s, n=nmo_not_converged, &
1336 6 : sym=dbcsr_type_no_symmetry)
1337 :
1338 6 : CALL dbcsr_init_p(matrix_sc)
1339 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_sc, template=matrix_s, n=nmo_not_converged, &
1340 6 : sym=dbcsr_type_no_symmetry)
1341 :
1342 6 : CALL dbcsr_init_p(matrix_z)
1343 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_z, template=matrix_s, n=nmo_not_converged, &
1344 6 : sym=dbcsr_type_no_symmetry)
1345 :
1346 6 : CALL copy_fm_to_dbcsr(mo_notconv_fm, mo_notconv) !fm->dbcsr
1347 :
1348 : CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, mo_notconv, 0.0_dp, matrix_hc, &
1349 6 : last_column=nmo_not_converged)
1350 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_notconv, 0.0_dp, matrix_sc, &
1351 6 : last_column=nmo_not_converged)
1352 :
1353 6 : CALL dbcsr_copy(matrix_z, matrix_sc)
1354 6 : CALL dbcsr_scale_by_vector(matrix_z, eig_not_conv, side='right')
1355 6 : CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
1356 :
1357 6 : DEALLOCATE (eig_not_conv)
1358 :
1359 : ! matrix_mm
1360 6 : CALL dbcsr_init_p(matrix_mm)
1361 : CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo_not_converged, n=nmo_not_converged, &
1362 6 : sym=dbcsr_type_no_symmetry)
1363 :
1364 : CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_notconv, matrix_hc, 0.0_dp, matrix_mm, &
1365 18 : last_column=nmo_not_converged)
1366 :
1367 : ELSE
1368 58 : mo_notconv => mo_coeff_b
1369 58 : mo_notconv_fm => mo_coeff
1370 58 : c_out => matrix_h
1371 : END IF
1372 :
1373 : ! allocate matrix_pz using as template matrix_z
1374 64 : CALL dbcsr_init_p(matrix_pz)
1375 : CALL dbcsr_create(matrix_pz, template=matrix_z, &
1376 : name="matrix_pz", &
1377 64 : matrix_type=dbcsr_type_no_symmetry)
1378 :
1379 64 : IF (do_apply_preconditioner) THEN
1380 56 : IF (preconditioner%in_use /= 0) THEN
1381 56 : CALL apply_preconditioner(preconditioner, matrix_z, matrix_pz)
1382 : ELSE
1383 0 : CALL dbcsr_copy(matrix_pz, matrix_z)
1384 : END IF
1385 : ELSE
1386 8 : CALL dbcsr_copy(matrix_pz, matrix_z)
1387 : END IF
1388 :
1389 : !allocate NMOxNMO full matrices
1390 64 : nmat = nmo_not_converged
1391 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat, ncol_global=nmat, &
1392 : context=mo_coeff%matrix_struct%context, &
1393 64 : para_env=mo_coeff%matrix_struct%para_env)
1394 64 : CALL cp_fm_create(matrix_mm_fm, fm_struct_tmp, name="m_tmp_mxm")
1395 64 : CALL cp_fm_create(matrix_mmt_fm, fm_struct_tmp, name="mt_tmp_mxm")
1396 64 : CALL cp_fm_struct_release(fm_struct_tmp)
1397 :
1398 : !allocate 2NMOx2NMO full matrices
1399 64 : nmat2 = 2*nmo_not_converged
1400 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
1401 : context=mo_coeff%matrix_struct%context, &
1402 64 : para_env=mo_coeff%matrix_struct%para_env)
1403 :
1404 64 : CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
1405 64 : CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
1406 64 : CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
1407 64 : CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
1408 192 : ALLOCATE (evals(nmat2))
1409 64 : CALL cp_fm_struct_release(fm_struct_tmp)
1410 :
1411 : ! compute CSC
1412 64 : CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
1413 : ! compute CHC
1414 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1415 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1, 1)
1416 :
1417 : ! compute the bottom left ZSC (top right is transpose)
1418 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=nmat)
1419 : ! set the bottom left part of S[C,Z] block matrix ZSC
1420 : !copy sparse to full
1421 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1422 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
1423 64 : CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
1424 64 : CALL cp_fm_to_fm_submat(matrix_mmt_fm, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
1425 :
1426 : ! compute the bottom left ZHC (top right is transpose)
1427 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=nmat)
1428 : ! set the bottom left part of S[C,Z] block matrix ZHC
1429 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1430 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
1431 64 : CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
1432 64 : CALL cp_fm_to_fm_submat(matrix_mmt_fm, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
1433 :
1434 64 : CALL cp_fm_release(matrix_mmt_fm)
1435 :
1436 : ! (reuse matrix_sc and matrix_hc to computr HZ and SZ)
1437 64 : CALL dbcsr_get_info(matrix_pz, nfullrows_total=n, nfullcols_total=k)
1438 64 : CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, matrix_pz, 0.0_dp, matrix_hc, last_column=k)
1439 64 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, matrix_pz, 0.0_dp, matrix_sc, last_column=k)
1440 :
1441 : ! compute the bottom right ZSZ
1442 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=k)
1443 : ! set the bottom right part of S[C,Z] block matrix ZSZ
1444 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1445 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
1446 :
1447 : ! compute the bottom right ZHZ
1448 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
1449 : ! set the bottom right part of H[C,Z] block matrix ZHZ
1450 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1451 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
1452 :
1453 64 : CALL dbcsr_release_p(matrix_mm)
1454 64 : CALL dbcsr_release_p(matrix_sc)
1455 64 : CALL dbcsr_release_p(matrix_hc)
1456 :
1457 64 : CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
1458 :
1459 : ! allocate two (nao x nmat) full matrix
1460 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmat, &
1461 : context=mo_coeff%matrix_struct%context, &
1462 64 : para_env=mo_coeff%matrix_struct%para_env)
1463 64 : CALL cp_fm_create(matrix_nm_fm, fm_struct_tmp, name="m_nxm")
1464 64 : CALL cp_fm_create(matrix_z_fm, fm_struct_tmp, name="m_nxm")
1465 64 : CALL cp_fm_struct_release(fm_struct_tmp)
1466 :
1467 64 : CALL copy_dbcsr_to_fm(matrix_pz, matrix_z_fm)
1468 : ! extract egenvectors
1469 64 : CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1, 1, 1, 1)
1470 64 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, mo_notconv_fm, matrix_mm_fm, 0.0_dp, matrix_nm_fm)
1471 64 : CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1 + nmat, 1, 1, 1)
1472 64 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, matrix_z_fm, matrix_mm_fm, 1.0_dp, matrix_nm_fm)
1473 :
1474 64 : CALL dbcsr_release_p(matrix_z)
1475 64 : CALL dbcsr_release_p(matrix_pz)
1476 64 : CALL cp_fm_release(matrix_z_fm)
1477 64 : CALL cp_fm_release(s_block)
1478 64 : CALL cp_fm_release(h_block)
1479 64 : CALL cp_fm_release(w_block)
1480 64 : CALL cp_fm_release(v_block)
1481 64 : CALL cp_fm_release(matrix_mm_fm)
1482 :
1483 : ! in case some vector are already converged only a subset of vectors are copied in the MOS
1484 64 : IF (nmo_converged > 0) THEN
1485 6 : jj = 1
1486 24 : DO j = 1, nset_not_conv
1487 18 : i_first = inotconv_set(j, 1)
1488 18 : i_last = inotconv_set(j, 2)
1489 18 : n = i_last - i_first + 1
1490 18 : CALL cp_fm_to_fm_submat(matrix_nm_fm, mo_coeff, nao, n, 1, jj, 1, i_first)
1491 186 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
1492 24 : jj = jj + n
1493 : END DO
1494 6 : DEALLOCATE (iconv_set)
1495 6 : DEALLOCATE (inotconv_set)
1496 :
1497 6 : CALL dbcsr_release_p(mo_notconv)
1498 6 : CALL dbcsr_release_p(c_out)
1499 6 : CALL cp_fm_release(mo_notconv_fm)
1500 6 : DEALLOCATE (mo_notconv_fm)
1501 : ELSE
1502 58 : CALL cp_fm_to_fm(matrix_nm_fm, mo_coeff)
1503 2146 : eigenvalues(1:nmo) = evals(1:nmo)
1504 : END IF
1505 64 : DEALLOCATE (evals)
1506 :
1507 64 : CALL cp_fm_release(matrix_nm_fm)
1508 64 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
1509 :
1510 64 : t2 = m_walltime()
1511 64 : IF (output_unit > 0) THEN
1512 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1513 0 : iteration, nmo_converged, max_norm, min_norm, t2 - t1
1514 : END IF
1515 536 : t1 = m_walltime()
1516 :
1517 : END DO ! iteration
1518 :
1519 64 : DEALLOCATE (ritz_coeff)
1520 64 : DEALLOCATE (iconv)
1521 64 : DEALLOCATE (inotconv)
1522 64 : DEALLOCATE (vnorm)
1523 :
1524 64 : CALL timestop(handle)
1525 :
1526 192 : END SUBROUTINE generate_extended_space_sparse
1527 :
1528 : ! **************************************************************************************************
1529 :
1530 : ! **************************************************************************************************
1531 : !> \brief ...
1532 : !> \param s_block ...
1533 : !> \param h_block ...
1534 : !> \param v_block ...
1535 : !> \param w_block ...
1536 : !> \param evals ...
1537 : !> \param ndim ...
1538 : ! **************************************************************************************************
1539 148 : SUBROUTINE reduce_extended_space(s_block, h_block, v_block, w_block, evals, ndim)
1540 :
1541 : TYPE(cp_fm_type), INTENT(IN) :: s_block, h_block, v_block, w_block
1542 : REAL(dp), DIMENSION(:) :: evals
1543 : INTEGER :: ndim
1544 :
1545 : CHARACTER(len=*), PARAMETER :: routineN = 'reduce_extended_space'
1546 :
1547 : INTEGER :: handle, info
1548 :
1549 148 : CALL timeset(routineN, handle)
1550 :
1551 148 : CALL cp_fm_to_fm(s_block, w_block)
1552 148 : CALL cp_fm_cholesky_decompose(s_block, info_out=info)
1553 148 : IF (info == 0) THEN
1554 148 : CALL cp_fm_triangular_invert(s_block)
1555 148 : CALL cp_fm_cholesky_restore(H_block, ndim, S_block, w_block, "MULTIPLY", pos="RIGHT")
1556 148 : CALL cp_fm_cholesky_restore(w_block, ndim, S_block, H_block, "MULTIPLY", pos="LEFT", transa="T")
1557 148 : CALL choose_eigv_solver(H_block, w_block, evals)
1558 148 : CALL cp_fm_cholesky_restore(w_block, ndim, S_block, v_block, "MULTIPLY")
1559 : ELSE
1560 : ! S^(-1/2)
1561 0 : CALL cp_fm_power(w_block, s_block, -0.5_dp, 1.0E-5_dp, info)
1562 0 : CALL cp_fm_to_fm(w_block, s_block)
1563 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, H_block, s_block, 0.0_dp, w_block)
1564 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, H_block)
1565 0 : CALL choose_eigv_solver(H_block, w_block, evals)
1566 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, v_block)
1567 : END IF
1568 :
1569 148 : CALL timestop(handle)
1570 :
1571 148 : END SUBROUTINE reduce_extended_space
1572 :
1573 : END MODULE qs_scf_block_davidson
|
