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