(git:f2099e5)
Loading...
Searching...
No Matches
qs_scf_block_davidson.F
Go to the documentation of this file.
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! **************************************************************************************************
18MODULE 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
85CONTAINS
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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
110 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 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 CALL timeset(routinen, handle)
123
124 NULLIFY (mo_coeff, mo_coeff_b, eigenvalues)
125
126 do_apply_preconditioner = .false.
127 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 nao=nao, nmo=nmo, homo=homo)
130 IF (do_apply_preconditioner) THEN
131 max_iter = bdav_env%max_iter
132 ELSE
133 max_iter = 1
134 END IF
135
136 NULLIFY (c_z, c_pz)
137 NULLIFY (evals, eig_not_conv)
138 t1 = m_walltime()
139 IF (output_unit > 0) THEN
140 WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
141 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
142 END IF
143
144 ALLOCATE (iconv(nmo))
145 ALLOCATE (inotconv(nmo))
146 ALLOCATE (ritz_coeff(nmo))
147 ALLOCATE (vnorm(nmo))
148
149 converged = .false.
150 DO iter = 1, max_iter
151
152 ! compute Ritz values
153 ritz_coeff = 0.0_dp
154 CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc")
155 CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo)
156 CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc")
157 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 para_env=mo_coeff%matrix_struct%para_env)
162 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
163 CALL cp_fm_struct_release(fm_struct_tmp)
164
165 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp)
166 CALL cp_fm_get_diag(m_tmp, ritz_coeff)
167 CALL cp_fm_release(m_tmp)
168
169 ! Check for converged eigenvectors
170 c_z => bdav_env%matrix_z
171 c_pz => bdav_env%matrix_pz
172 CALL cp_fm_to_fm(m_sc, c_z)
173 CALL cp_fm_column_scale(c_z, ritz_coeff)
174 CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc)
175 CALL cp_fm_vectorsnorm(c_z, vnorm)
176
177 nmo_converged = 0
178 nmo_not_converged = 0
179 max_norm = 0.0_dp
180 min_norm = 1.e10_dp
181 DO imo = 1, nmo
182 max_norm = max(max_norm, vnorm(imo))
183 min_norm = min(min_norm, vnorm(imo))
184 END DO
185 iconv = 0
186 inotconv = 0
187 DO imo = 1, nmo
188 IF (vnorm(imo) <= bdav_env%eps_iter) THEN
189 nmo_converged = nmo_converged + 1
190 iconv(nmo_converged) = imo
191 ELSE
192 nmo_not_converged = nmo_not_converged + 1
193 inotconv(nmo_not_converged) = imo
194 END IF
195 END DO
196
197 IF (nmo_converged > 0) THEN
198 ALLOCATE (iconv_set(nmo_converged, 2))
199 ALLOCATE (inotconv_set(nmo_not_converged, 2))
200 i_last = iconv(1)
201 nset = 0
202 DO j = 1, nmo_converged
203 imo = iconv(j)
204
205 IF (imo == i_last + 1) THEN
206 i_last = imo
207 iconv_set(nset, 2) = imo
208 ELSE
209 i_last = imo
210 nset = nset + 1
211 iconv_set(nset, 1) = imo
212 iconv_set(nset, 2) = imo
213 END IF
214 END DO
215 nset_conv = nset
216
217 i_last = inotconv(1)
218 nset = 0
219 DO j = 1, nmo_not_converged
220 imo = inotconv(j)
221
222 IF (imo == i_last + 1) THEN
223 i_last = imo
224 inotconv_set(nset, 2) = imo
225 ELSE
226 i_last = imo
227 nset = nset + 1
228 inotconv_set(nset, 1) = imo
229 inotconv_set(nset, 2) = imo
230 END IF
231 END DO
232 nset_not_conv = nset
233 CALL cp_fm_release(m_sc)
234 CALL cp_fm_release(m_hc)
235 NULLIFY (c_z, c_pz)
236 END IF
237
238 IF (real(nmo_converged, dp)/real(nmo, dp) > bdav_env%conv_percent) THEN
239 converged = .true.
240 DEALLOCATE (iconv_set)
241 DEALLOCATE (inotconv_set)
242 t2 = m_walltime()
243 IF (output_unit > 0) THEN
244 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
245 iter, nmo_converged, max_norm, min_norm, t2 - t1
246
247 WRITE (output_unit, *) " Reached convergence in ", iter, &
248 " Davidson iterations"
249 END IF
250
251 EXIT
252 END IF
253
254 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 para_env=mo_coeff%matrix_struct%para_env)
258 !allocate h_fm
259 CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp")
260 !allocate s_fm
261 CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp")
262 !copy matrix_h in h_fm
263 CALL copy_dbcsr_to_fm(matrix_h, h_fm)
264 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 CALL copy_dbcsr_to_fm(matrix_s, s_fm)
269
270 !allocate c_out
271 CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp")
272 ! set c_out to zero
273 CALL cp_fm_set_all(c_out, 0.0_dp)
274 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 para_env=mo_coeff%matrix_struct%para_env)
280 CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv")
281 CALL cp_fm_set_all(c_conv, 0.0_dp)
282 !allocate m_tmp
283 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc")
284 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 para_env=mo_coeff%matrix_struct%para_env)
291 CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv")
292 CALL cp_fm_set_all(c_notconv, 0.0_dp)
293 IF (nmo_converged > 0) THEN
294 CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc")
295 CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc")
296 !allocate c_z
297 ALLOCATE (c_z, c_pz)
298 CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z")
299 CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz")
300 CALL cp_fm_set_all(c_z, 0.0_dp)
301
302 ! sum contributions to c_out
303 jj = 1
304 DO j = 1, nset_conv
305 i_first = iconv_set(j, 1)
306 i_last = iconv_set(j, 2)
307 n = i_last - i_first + 1
308 CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj)
309 jj = jj + n
310 END DO
311 CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp)
312 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 lambda = 100.0_dp*abs(eigenvalues(homo))
316 CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm)
317 CALL cp_fm_release(m_tmp)
318 CALL cp_fm_release(h_fm)
319
320 END IF
321
322 !allocate m_tmp
323 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm")
324 CALL cp_fm_struct_release(fm_struct_tmp)
325 IF (nmo_converged > 0) THEN
326 ALLOCATE (eig_not_conv(nmo_not_converged))
327 jj = 1
328 DO j = 1, nset_not_conv
329 i_first = inotconv_set(j, 1)
330 i_last = inotconv_set(j, 2)
331 n = i_last - i_first + 1
332 CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj)
333 eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
334 jj = jj + n
335 END DO
336 CALL parallel_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc)
337 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 CALL cp_fm_to_fm(m_sc, m_tmp)
340 CALL cp_fm_column_scale(m_tmp, eig_not_conv)
341 CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc)
342 DEALLOCATE (eig_not_conv)
343 CALL cp_fm_to_fm(m_tmp, c_z)
344 ELSE
345 CALL cp_fm_to_fm(mo_coeff, c_notconv)
346 END IF
347
348 !preconditioner
349 IF (do_apply_preconditioner) THEN
350 IF (preconditioner%in_use /= 0) THEN
351 CALL apply_preconditioner(preconditioner, c_z, c_pz)
352 ELSE
353 CALL cp_fm_to_fm(c_z, c_pz)
354 END IF
355 ELSE
356 CALL cp_fm_to_fm(c_z, c_pz)
357 END IF
358 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 para_env=mo_coeff%matrix_struct%para_env)
363
364 CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm")
365 CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm")
366 CALL cp_fm_struct_release(fm_struct_tmp)
367
368 nmat = nmo_not_converged
369 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 para_env=mo_coeff%matrix_struct%para_env)
373
374 CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
375 CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
376 CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
377 CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
378 ALLOCATE (evals(nmat2))
379
380 CALL cp_fm_struct_release(fm_struct_tmp)
381
382 ! compute CSC
383 CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
384
385 ! compute CHC
386 CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp)
387 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1)
388
389 ! compute ZSC
390 CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
391 CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
392 CALL cp_fm_transpose(m_tmp, mt_tmp)
393 CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
394 ! compute ZHC
395 CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
396 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
397 CALL cp_fm_transpose(m_tmp, mt_tmp)
398 CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
399
400 CALL cp_fm_release(mt_tmp)
401
402 ! reuse m_sc and m_hc to computr HZ and SZ
403 IF (nmo_converged > 0) THEN
404 CALL parallel_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc)
405 CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc)
406
407 CALL cp_fm_release(c_out)
408 CALL cp_fm_release(c_conv)
409 CALL cp_fm_release(s_fm)
410 ELSE
411 CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo)
412 CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo)
413 END IF
414
415 ! compute ZSZ
416 CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
417 CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
418 ! compute ZHZ
419 CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
420 CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
421
422 CALL cp_fm_release(m_sc)
423
424 ! solution of the reduced eigenvalues problem
425 CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
426
427 ! extract egenvectors
428 CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1)
429 CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc)
430 CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1)
431 CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc)
432
433 CALL cp_fm_release(m_tmp)
434
435 CALL cp_fm_release(c_notconv)
436 CALL cp_fm_release(s_block)
437 CALL cp_fm_release(h_block)
438 CALL cp_fm_release(w_block)
439 CALL cp_fm_release(v_block)
440
441 IF (nmo_converged > 0) THEN
442 CALL cp_fm_release(c_z)
443 CALL cp_fm_release(c_pz)
444 DEALLOCATE (c_z, c_pz)
445 jj = 1
446 DO j = 1, nset_not_conv
447 i_first = inotconv_set(j, 1)
448 i_last = inotconv_set(j, 2)
449 n = i_last - i_first + 1
450 CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first)
451 eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
452 jj = jj + n
453 END DO
454 DEALLOCATE (iconv_set)
455 DEALLOCATE (inotconv_set)
456 ELSE
457 CALL cp_fm_to_fm(m_hc, mo_coeff)
458 eigenvalues(1:nmo) = evals(1:nmo)
459 END IF
460 DEALLOCATE (evals)
461
462 CALL cp_fm_release(m_hc)
463
464 CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
465
466 t2 = m_walltime()
467 IF (output_unit > 0) THEN
468 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
469 iter, nmo_converged, max_norm, min_norm, t2 - t1
470 END IF
471 t1 = m_walltime()
472
473 END DO ! iter
474
475 DEALLOCATE (iconv)
476 DEALLOCATE (inotconv)
477 DEALLOCATE (ritz_coeff)
478 DEALLOCATE (vnorm)
479
480 CALL timestop(handle)
481 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdiag, scaling
515 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
519 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm, vnorm_pz
525 REAL(kind=dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals, &
526 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 CALL timeset(routinen, handle)
544
545 NULLIFY (eigenvalues, evals, eig_not_conv, cmo, occupation, fm_struct_tmp)
546
547 do_apply_preconditioner = .false.
548 IF (PRESENT(preconditioner)) do_apply_preconditioner = .true.
549 CALL get_mo_set(mo_set=mos, cmo_coeff=cmo, eigenvalues=eigenvalues, &
550 occupation_numbers=occupation, nao=nao, nmo=nmo, homo=homo)
551 IF (do_apply_preconditioner) THEN
552 max_iter = bdav_env%max_iter
553 ELSE
554 max_iter = 1
555 END IF
556
557 t1 = m_walltime()
558 IF (output_unit > 0) THEN
559 WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
560 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
561 END IF
562
563 ALLOCATE (iconv(nmo))
564 ALLOCATE (inotconv(nmo))
565 ALLOCATE (ritz_coeff(nmo))
566 ALLOCATE (vnorm(nmo))
567 ALLOCATE (cdiag(nmo))
568
569 converged = .false.
570 exhausted = .false.
571 DO iter = 1, max_iter
572
573 ! compute Ritz values
574 ritz_coeff = 0.0_dp
575 CALL cp_cfm_create(m_hc, cmo%matrix_struct, name="hc")
576 CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, matrix_h, cmo, z_zero, m_hc)
577 CALL cp_cfm_create(m_sc, cmo%matrix_struct, name="sc")
578 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 para_env=cmo%matrix_struct%para_env)
583 CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
584 CALL cp_fm_struct_release(fm_struct_tmp)
585
586 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmo, m_hc, z_zero, m_tmp)
587 CALL cp_cfm_get_diag(m_tmp, cdiag)
588 ritz_coeff(1:nmo) = real(cdiag(1:nmo), kind=dp)
589 CALL cp_cfm_release(m_tmp)
590
591 ! Check for converged eigenvectors: residual C_z = S*C*diag(ritz) - H*C
592 CALL cp_cfm_create(c_z, cmo%matrix_struct, name="z")
593 CALL cp_cfm_to_cfm(m_sc, c_z)
594 IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
595 ALLOCATE (scaling(nmo))
596 scaling(:) = cmplx(ritz_coeff, 0.0_dp, kind=dp)
597 CALL cp_cfm_column_scale(c_z, scaling)
598 CALL cp_cfm_scale_and_add(-z_one, c_z, z_one, m_hc)
599 CALL cp_cfm_vectorsnorm(c_z, vnorm)
600
601 nmo_converged = 0
602 nmo_not_converged = 0
603 max_norm = 0.0_dp
604 min_norm = 1.e10_dp
605 DO imo = 1, nmo
606 max_norm = max(max_norm, vnorm(imo))
607 min_norm = min(min_norm, vnorm(imo))
608 END DO
609 iconv = 0
610 inotconv = 0
611 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 IF (occupation(imo) < occ_tol) THEN
617 eps_iter_col = eps_iter_empty
618 ELSE
619 eps_iter_col = eps_iter
620 END IF
621 IF (vnorm(imo) <= eps_iter_col) THEN
622 nmo_converged = nmo_converged + 1
623 iconv(nmo_converged) = imo
624 ELSE
625 nmo_not_converged = nmo_not_converged + 1
626 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 IF (iter == 1 .AND. nmo_not_converged == 0) THEN
641 nmo_not_converged = nmo
642 nmo_converged = 0
643 iconv = 0
644 DO imo = 1, nmo
645 inotconv(imo) = imo
646 END DO
647 END IF
648
649 IF (nmo_converged > 0) THEN
650 ALLOCATE (iconv_set(nmo_converged, 2))
651 ALLOCATE (inotconv_set(nmo_not_converged, 2))
652 i_last = iconv(1)
653 nset = 0
654 DO j = 1, nmo_converged
655 imo = iconv(j)
656
657 IF (imo == i_last + 1) THEN
658 i_last = imo
659 iconv_set(nset, 2) = imo
660 ELSE
661 i_last = imo
662 nset = nset + 1
663 iconv_set(nset, 1) = imo
664 iconv_set(nset, 2) = imo
665 END IF
666 END DO
667 nset_conv = nset
668
669 i_last = inotconv(1)
670 nset = 0
671 DO j = 1, nmo_not_converged
672 imo = inotconv(j)
673
674 IF (imo == i_last + 1) THEN
675 i_last = imo
676 inotconv_set(nset, 2) = imo
677 ELSE
678 i_last = imo
679 nset = nset + 1
680 inotconv_set(nset, 1) = imo
681 inotconv_set(nset, 2) = imo
682 END IF
683 END DO
684 nset_not_conv = nset
685 CALL cp_cfm_release(m_sc)
686 CALL cp_cfm_release(m_hc)
687 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 IF (iter > 1 .AND. real(nmo_converged, dp)/real(nmo, dp) > bdav_env%conv_percent) THEN
699 converged = .true.
700 DEALLOCATE (iconv_set)
701 DEALLOCATE (inotconv_set)
702 t2 = m_walltime()
703 IF (output_unit > 0) THEN
704 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
705 iter, nmo_converged, max_norm, min_norm, t2 - t1
706
707 WRITE (output_unit, *) " Reached convergence in ", iter, &
708 " Davidson iterations"
709 END IF
710
711 EXIT
712 END IF
713
714 ncol_z = nmo_not_converged
715
716 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 para_env=cmo%matrix_struct%para_env)
721 CALL cp_cfm_create(h_fm, fm_struct_tmp, name="hf")
722 CALL cp_cfm_create(c_out, fm_struct_tmp, name="cout")
723 CALL cp_fm_struct_release(fm_struct_tmp)
724 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 para_env=cmo%matrix_struct%para_env)
730 CALL cp_cfm_create(c_conv, fm_struct_tmp, name="cconv")
731 CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="scconv")
732 CALL cp_fm_struct_release(fm_struct_tmp)
733 jj = 1
734 DO j = 1, nset_conv
735 i_first = iconv_set(j, 1)
736 i_last = iconv_set(j, 2)
737 n = i_last - i_first + 1
738 ALLOCATE (submat_buffer(nao, n))
739 CALL cp_cfm_get_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first, &
740 n_rows=nao, n_cols=n)
741 CALL cp_cfm_set_submatrix(c_conv, submat_buffer, start_row=1, start_col=jj)
742 DEALLOCATE (submat_buffer)
743 jj = jj + n
744 END DO
745 CALL cp_cfm_gemm('N', 'N', nao, nmo_converged, nao, z_one, matrix_s, c_conv, &
746 z_zero, m_tmp)
747 CALL cp_cfm_gemm('N', 'C', nao, nao, nmo_converged, z_one, m_tmp, m_tmp, z_zero, c_out)
748 CALL cp_cfm_release(m_tmp)
749 CALL cp_cfm_release(c_conv)
750
751 ! project c_out out of H
752 lambda = 100.0_dp*abs(eigenvalues(homo))
753 lambda_c = cmplx(lambda, 0.0_dp, kind=dp)
754 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 para_env=cmo%matrix_struct%para_env)
761 CALL cp_cfm_create(c_notconv, fm_struct_tmp, name="c_notconv")
762 CALL cp_fm_struct_release(fm_struct_tmp)
763
764 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 para_env=cmo%matrix_struct%para_env)
769 CALL cp_cfm_create(m_hc, fm_struct_tmp, name="m_hc")
770 CALL cp_cfm_create(m_sc, fm_struct_tmp, name="m_sc")
771 CALL cp_cfm_create(c_z, fm_struct_tmp, name="c_z")
772 CALL cp_fm_struct_release(fm_struct_tmp)
773
774 ALLOCATE (eig_not_conv(ncol_z))
775 jj = 1
776 DO j = 1, nset_not_conv
777 i_first = inotconv_set(j, 1)
778 i_last = inotconv_set(j, 2)
779 n = i_last - i_first + 1
780 ALLOCATE (submat_buffer(nao, n))
781 CALL cp_cfm_get_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first, &
782 n_rows=nao, n_cols=n)
783 CALL cp_cfm_set_submatrix(c_notconv, submat_buffer, start_row=1, start_col=jj)
784 DEALLOCATE (submat_buffer)
785 eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
786 jj = jj + n
787 END DO
788 ! extend the subspace using only the not converged vectors
789 CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, h_fm, c_notconv, z_zero, m_hc)
790 CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_s, c_notconv, z_zero, m_sc)
791 CALL cp_cfm_to_cfm(m_sc, c_z)
792 IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
793 ALLOCATE (scaling(ncol_z))
794 scaling(:) = cmplx(eig_not_conv, 0.0_dp, kind=dp)
795 CALL cp_cfm_column_scale(c_z, scaling)
796 CALL cp_cfm_scale_and_add(-z_one, c_z, z_one, m_hc)
797 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 CALL cp_cfm_to_cfm(cmo, c_notconv)
803 END IF
804
805 ! preconditioner
806 CALL cp_cfm_create(c_mz, c_z%matrix_struct, name="pz")
807 IF (do_apply_preconditioner) THEN
808 IF (preconditioner%in_use /= 0) THEN
809 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 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 CALL cp_cfm_create(c_pzw, cmo%matrix_struct, name="pzw")
823 CALL cp_cfm_create(c_pzo, cmo%matrix_struct, name="pzo")
824 CALL cp_cfm_set_all(c_pzw, z_zero, z_zero)
825 ALLOCATE (colbuf(nao, ncol_z))
826 jj = 1
827 DO j = 1, nset_not_conv
828 i_first = inotconv_set(j, 1)
829 i_last = inotconv_set(j, 2)
830 n = i_last - i_first + 1
831 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 start_col=i_first)
834 jj = jj + n
835 END DO
836 CALL apply_preconditioner(preconditioner, c_pzw, c_pzo)
837 jj = 1
838 DO j = 1, nset_not_conv
839 i_first = inotconv_set(j, 1)
840 i_last = inotconv_set(j, 2)
841 n = i_last - i_first + 1
842 CALL cp_cfm_get_submatrix(c_pzo, colbuf(:, 1:n), n_rows=nao, n_cols=n, &
843 start_col=i_first)
844 CALL cp_cfm_set_submatrix(c_mz, colbuf(:, 1:n), n_rows=nao, n_cols=n, start_col=jj)
845 jj = jj + n
846 END DO
847 DEALLOCATE (colbuf)
848 CALL cp_cfm_release(c_pzw)
849 CALL cp_cfm_release(c_pzo)
850 END IF
851 ELSE
852 CALL cp_cfm_to_cfm(c_z, c_mz)
853 END IF
854 ELSE
855 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 ALLOCATE (vnorm_pz(ncol_z))
871 CALL cp_cfm_vectorsnorm(c_mz, vnorm_pz)
872 IF (any(vnorm_pz < eps_exhaust)) THEN
873 exhausted = .true.
874 ELSE
875 IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
876 ALLOCATE (scaling(ncol_z))
877 scaling(:) = cmplx(1.0_dp/vnorm_pz, 0.0_dp, kind=dp)
878 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 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 para_env=cmo%matrix_struct%para_env)
885 CALL cp_cfm_create(m_tmp, fm_struct_tmp, name="proj_coef")
886 CALL cp_fm_struct_release(fm_struct_tmp)
887 CALL cp_cfm_gemm('C', 'N', nmo, ncol_z, nao, z_one, cmo, c_z, z_zero, m_tmp)
888 CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nmo, -z_one, cmo, m_tmp, z_one, c_mz)
889 CALL cp_cfm_release(m_tmp)
890
891 CALL cp_cfm_vectorsnorm(c_mz, vnorm_pz)
892 IF (any(vnorm_pz < eps_exhaust)) THEN
893 exhausted = .true.
894 ELSE
895 IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
896 ALLOCATE (scaling(ncol_z))
897 scaling(:) = cmplx(1.0_dp/vnorm_pz, 0.0_dp, kind=dp)
898 CALL cp_cfm_column_scale(c_mz, scaling)
899 END IF
900 END IF
901 DEALLOCATE (vnorm_pz)
902 IF (exhausted) EXIT
903
904 nmat = ncol_z
905 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 para_env=cmo%matrix_struct%para_env)
909
910 CALL cp_cfm_create(s_block, fm_struct_tmp, name="sb")
911 CALL cp_cfm_create(h_block, fm_struct_tmp, name="hb")
912 CALL cp_cfm_create(v_block, fm_struct_tmp, name="vb")
913 CALL cp_cfm_create(w_block, fm_struct_tmp, name="wb")
914 ALLOCATE (evals(nmat2))
915 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 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 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 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 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 IF (nmo_converged > 0) THEN
935 CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, h_fm, c_mz, z_zero, m_hc)
936 CALL cp_cfm_release(c_out)
937 CALL cp_cfm_release(h_fm)
938 ELSE
939 CALL cp_cfm_gemm('N', 'N', nao, ncol_z, nao, z_one, matrix_h, c_mz, z_zero, m_hc)
940 END IF
941 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 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 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 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 c_first_row=1 + nmat, c_first_col=1 + nmat)
960
961 CALL cp_cfm_release(m_sc)
962
963 ! solution of the reduced generalized eigenproblem
964 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 para_env=cmo%matrix_struct%para_env)
970 CALL cp_cfm_create(c_new, fm_struct_tmp, name="c_new")
971 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 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 b_first_row=1 + nmat, b_first_col=1)
976
977 CALL cp_cfm_release(m_hc)
978 CALL cp_cfm_release(c_mz)
979 CALL cp_cfm_release(c_z)
980 CALL cp_cfm_release(c_notconv)
981 CALL cp_cfm_release(s_block)
982 CALL cp_cfm_release(h_block)
983 CALL cp_cfm_release(w_block)
984 CALL cp_cfm_release(v_block)
985
986 IF (nmo_converged > 0) THEN
987 jj = 1
988 DO j = 1, nset_not_conv
989 i_first = inotconv_set(j, 1)
990 i_last = inotconv_set(j, 2)
991 n = i_last - i_first + 1
992 ALLOCATE (submat_buffer(nao, n))
993 CALL cp_cfm_get_submatrix(c_new, submat_buffer, start_row=1, start_col=jj, &
994 n_rows=nao, n_cols=n)
995 CALL cp_cfm_set_submatrix(cmo, submat_buffer, start_row=1, start_col=i_first)
996 DEALLOCATE (submat_buffer)
997 eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
998 jj = jj + n
999 END DO
1000 DEALLOCATE (iconv_set)
1001 DEALLOCATE (inotconv_set)
1002 ELSE
1003 CALL cp_cfm_to_cfm(c_new, cmo)
1004 eigenvalues(1:nmo) = evals(1:nmo)
1005 END IF
1006 DEALLOCATE (evals)
1007 CALL cp_cfm_release(c_new)
1008
1009 t2 = m_walltime()
1010 IF (output_unit > 0) THEN
1011 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1012 iter, nmo_converged, max_norm, min_norm, t2 - t1
1013 END IF
1014 t1 = m_walltime()
1015
1016 END DO ! iter
1017
1018 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 IF (nmo_converged > 0) THEN
1023 CALL cp_cfm_release(h_fm)
1024 CALL cp_cfm_release(c_out)
1025 END IF
1026 CALL cp_cfm_release(m_hc)
1027 CALL cp_cfm_release(m_sc)
1028 CALL cp_cfm_release(c_z)
1029 CALL cp_cfm_release(c_mz)
1030 CALL cp_cfm_release(c_notconv)
1031 IF (ALLOCATED(iconv_set)) DEALLOCATE (iconv_set)
1032 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 CALL cp_cfm_create(h_fm, matrix_h%matrix_struct, name="fb_h")
1038 CALL cp_cfm_create(c_out, matrix_h%matrix_struct, name="fb_s")
1039 CALL cp_cfm_create(c_new, matrix_h%matrix_struct, name="fb_work")
1040 CALL cp_cfm_to_cfm(matrix_h, h_fm)
1041 CALL cp_cfm_to_cfm(matrix_s, c_out)
1042 CALL cp_cfm_geeig(h_fm, c_out, cmo, eigenvalues, c_new)
1043 CALL cp_cfm_release(h_fm)
1044 CALL cp_cfm_release(c_out)
1045 CALL cp_cfm_release(c_new)
1046 END IF
1047
1048 DEALLOCATE (iconv)
1049 DEALLOCATE (inotconv)
1050 DEALLOCATE (ritz_coeff)
1051 DEALLOCATE (vnorm)
1052 DEALLOCATE (cdiag)
1053 IF (ALLOCATED(scaling)) DEALLOCATE (scaling)
1054
1055 CALL timestop(handle)
1056 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
1081 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
1082 LOGICAL :: converged, do_apply_preconditioner
1083 REAL(dp) :: lambda, max_norm, min_norm, t1, t2
1084 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, ritz_coeff, vnorm
1085 REAL(dp), DIMENSION(:), POINTER :: eigenvalues
1086 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 CALL timeset(routinen, handle)
1099
1100 do_apply_preconditioner = .false.
1101 IF (PRESENT(preconditioner)) do_apply_preconditioner = .true.
1102
1103 NULLIFY (mo_coeff, mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm)
1104 NULLIFY (mo_notconv_fm, mo_conv, mo_notconv, smo_conv, c_out)
1105 NULLIFY (fm_struct_tmp)
1106 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
1107 eigenvalues=eigenvalues, homo=homo, nao=nao, nmo=nmo)
1108 IF (do_apply_preconditioner) THEN
1109 max_iter = bdav_env%max_iter
1110 ELSE
1111 max_iter = 1
1112 END IF
1113
1114 t1 = m_walltime()
1115 IF (output_unit > 0) THEN
1116 WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
1117 " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
1118 END IF
1119
1120 ! Allocate array for Ritz values
1121 ALLOCATE (ritz_coeff(nmo))
1122 ALLOCATE (iconv(nmo))
1123 ALLOCATE (inotconv(nmo))
1124 ALLOCATE (vnorm(nmo))
1125
1126 converged = .false.
1127 DO iteration = 1, max_iter
1128 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 CALL dbcsr_init_p(matrix_hc)
1131 CALL dbcsr_create(matrix_hc, template=mo_coeff_b, &
1132 name="matrix_hc", &
1133 matrix_type=dbcsr_type_no_symmetry)
1134 CALL dbcsr_init_p(matrix_sc)
1135 CALL dbcsr_create(matrix_sc, template=mo_coeff_b, &
1136 name="matrix_sc", &
1137 matrix_type=dbcsr_type_no_symmetry)
1138
1139 CALL dbcsr_get_info(mo_coeff_b, nfullrows_total=n, nfullcols_total=k, group=group)
1140 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h, mo_coeff_b, 0.0_dp, matrix_hc, last_column=k)
1141 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 ritz_coeff = 0.0_dp
1145 ! Allocate Sparse matrices: nmoxnmo
1146 ! matrix_mm
1147
1148 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 sym=dbcsr_type_no_symmetry)
1151
1152 CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_coeff_b, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
1153 CALL dbcsr_get_diag(matrix_mm, ritz_coeff)
1154 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 CALL dbcsr_init_p(matrix_z)
1158 CALL dbcsr_create(matrix_z, template=mo_coeff_b, &
1159 name="matrix_z", &
1160 matrix_type=dbcsr_type_no_symmetry)
1161 CALL dbcsr_copy(matrix_z, matrix_sc)
1162 CALL dbcsr_scale_by_vector(matrix_z, ritz_coeff, side='right')
1163 CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
1164
1165 ! Compute the column norms of matrix_z.
1166 vnorm = 0.0_dp
1167 CALL dbcsr_iterator_start(iter, matrix_z)
1168 DO WHILE (dbcsr_iterator_blocks_left(iter))
1169 CALL dbcsr_iterator_next_block(iter, block=block, col_offset=col_offset)
1170 DO j = 1, SIZE(block, 2)
1171 vnorm(col_offset + j - 1) = vnorm(col_offset + j - 1) + sum(block(:, j)**2)
1172 END DO
1173 END DO
1174 CALL dbcsr_iterator_stop(iter)
1175 CALL group%sum(vnorm)
1176 vnorm = sqrt(vnorm)
1177
1178 ! Check for converged eigenvectors
1179 nmo_converged = 0
1180 nmo_not_converged = 0
1181 max_norm = 0.0_dp
1182 min_norm = 1.e10_dp
1183 DO imo = 1, nmo
1184 max_norm = max(max_norm, vnorm(imo))
1185 min_norm = min(min_norm, vnorm(imo))
1186 END DO
1187 iconv = 0
1188 inotconv = 0
1189
1190 DO imo = 1, nmo
1191 IF (vnorm(imo) <= bdav_env%eps_iter) THEN
1192 nmo_converged = nmo_converged + 1
1193 iconv(nmo_converged) = imo
1194 ELSE
1195 nmo_not_converged = nmo_not_converged + 1
1196 inotconv(nmo_not_converged) = imo
1197 END IF
1198 END DO
1199
1200 IF (nmo_converged > 0) THEN
1201 ALLOCATE (iconv_set(nmo_converged, 2))
1202 ALLOCATE (inotconv_set(nmo_not_converged, 2))
1203 i_last = iconv(1)
1204 nset = 0
1205 DO j = 1, nmo_converged
1206 imo = iconv(j)
1207
1208 IF (imo == i_last + 1) THEN
1209 i_last = imo
1210 iconv_set(nset, 2) = imo
1211 ELSE
1212 i_last = imo
1213 nset = nset + 1
1214 iconv_set(nset, 1) = imo
1215 iconv_set(nset, 2) = imo
1216 END IF
1217 END DO
1218 nset_conv = nset
1219
1220 i_last = inotconv(1)
1221 nset = 0
1222 DO j = 1, nmo_not_converged
1223 imo = inotconv(j)
1224
1225 IF (imo == i_last + 1) THEN
1226 i_last = imo
1227 inotconv_set(nset, 2) = imo
1228 ELSE
1229 i_last = imo
1230 nset = nset + 1
1231 inotconv_set(nset, 1) = imo
1232 inotconv_set(nset, 2) = imo
1233 END IF
1234 END DO
1235 nset_not_conv = nset
1236
1237 CALL dbcsr_release_p(matrix_hc)
1238 CALL dbcsr_release_p(matrix_sc)
1239 CALL dbcsr_release_p(matrix_z)
1240 CALL dbcsr_release_p(matrix_mm)
1241 END IF
1242
1243 IF (real(nmo_converged, dp)/real(nmo, dp) > bdav_env%conv_percent) THEN
1244 DEALLOCATE (iconv_set)
1245
1246 DEALLOCATE (inotconv_set)
1247
1248 converged = .true.
1249 t2 = m_walltime()
1250 IF (output_unit > 0) THEN
1251 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1252 iteration, nmo_converged, max_norm, min_norm, t2 - t1
1253
1254 WRITE (output_unit, *) " Reached convergence in ", iteration, &
1255 " Davidson iterations"
1256 END IF
1257
1258 EXIT
1259 END IF
1260
1261 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 para_env=mo_coeff%matrix_struct%para_env)
1267 CALL cp_fm_create(mo_conv_fm, fm_struct_tmp, name="mo_conv_fm")
1268
1269 CALL cp_fm_struct_release(fm_struct_tmp)
1270
1271 ! extract mo_conv from mo_coeff full matrix
1272 jj = 1
1273 DO j = 1, nset_conv
1274 i_first = iconv_set(j, 1)
1275 i_last = iconv_set(j, 2)
1276 n = i_last - i_first + 1
1277 CALL cp_fm_to_fm_submat(mo_coeff, mo_conv_fm, nao, n, 1, i_first, 1, jj)
1278 jj = jj + n
1279 END DO
1280
1281 ! allocate c_out sparse matrix, to project out the converged MOS
1282 CALL dbcsr_init_p(c_out)
1283 CALL dbcsr_create(c_out, template=matrix_s, &
1284 name="c_out", &
1285 matrix_type=dbcsr_type_symmetric)
1286
1287 ! allocate mo_conv sparse
1288 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 sym=dbcsr_type_no_symmetry)
1291
1292 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 sym=dbcsr_type_no_symmetry)
1295
1296 CALL copy_fm_to_dbcsr(mo_conv_fm, mo_conv) !fm->dbcsr
1297
1298 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_conv, 0.0_dp, smo_conv, last_column=nmo_converged)
1299 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 lambda = 100.0_dp*abs(eigenvalues(homo))
1302 CALL dbcsr_add(c_out, matrix_h, lambda, 1.0_dp)
1303
1304 CALL dbcsr_release_p(mo_conv)
1305 CALL dbcsr_release_p(smo_conv)
1306 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 para_env=mo_coeff%matrix_struct%para_env)
1312 ALLOCATE (mo_notconv_fm)
1313 CALL cp_fm_create(mo_notconv_fm, fm_struct_tmp, name="mo_notconv_fm")
1314 CALL cp_fm_struct_release(fm_struct_tmp)
1315
1316 ! extract mo_notconv from mo_coeff full matrix
1317 ALLOCATE (eig_not_conv(nmo_not_converged))
1318
1319 jj = 1
1320 DO j = 1, nset_not_conv
1321 i_first = inotconv_set(j, 1)
1322 i_last = inotconv_set(j, 2)
1323 n = i_last - i_first + 1
1324 CALL cp_fm_to_fm_submat(mo_coeff, mo_notconv_fm, nao, n, 1, i_first, 1, jj)
1325 eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
1326 jj = jj + n
1327 END DO
1328
1329 ! allocate mo_conv sparse
1330 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 sym=dbcsr_type_no_symmetry)
1333
1334 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 sym=dbcsr_type_no_symmetry)
1337
1338 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 sym=dbcsr_type_no_symmetry)
1341
1342 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 sym=dbcsr_type_no_symmetry)
1345
1346 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 last_column=nmo_not_converged)
1350 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_notconv, 0.0_dp, matrix_sc, &
1351 last_column=nmo_not_converged)
1352
1353 CALL dbcsr_copy(matrix_z, matrix_sc)
1354 CALL dbcsr_scale_by_vector(matrix_z, eig_not_conv, side='right')
1355 CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
1356
1357 DEALLOCATE (eig_not_conv)
1358
1359 ! matrix_mm
1360 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 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 last_column=nmo_not_converged)
1366
1367 ELSE
1368 mo_notconv => mo_coeff_b
1369 mo_notconv_fm => mo_coeff
1370 c_out => matrix_h
1371 END IF
1372
1373 ! allocate matrix_pz using as template matrix_z
1374 CALL dbcsr_init_p(matrix_pz)
1375 CALL dbcsr_create(matrix_pz, template=matrix_z, &
1376 name="matrix_pz", &
1377 matrix_type=dbcsr_type_no_symmetry)
1378
1379 IF (do_apply_preconditioner) THEN
1380 IF (preconditioner%in_use /= 0) THEN
1381 CALL apply_preconditioner(preconditioner, matrix_z, matrix_pz)
1382 ELSE
1383 CALL dbcsr_copy(matrix_pz, matrix_z)
1384 END IF
1385 ELSE
1386 CALL dbcsr_copy(matrix_pz, matrix_z)
1387 END IF
1388
1389 !allocate NMOxNMO full matrices
1390 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 para_env=mo_coeff%matrix_struct%para_env)
1394 CALL cp_fm_create(matrix_mm_fm, fm_struct_tmp, name="m_tmp_mxm")
1395 CALL cp_fm_create(matrix_mmt_fm, fm_struct_tmp, name="mt_tmp_mxm")
1396 CALL cp_fm_struct_release(fm_struct_tmp)
1397
1398 !allocate 2NMOx2NMO full matrices
1399 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 para_env=mo_coeff%matrix_struct%para_env)
1403
1404 CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
1405 CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
1406 CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
1407 CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
1408 ALLOCATE (evals(nmat2))
1409 CALL cp_fm_struct_release(fm_struct_tmp)
1410
1411 ! compute CSC
1412 CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
1413 ! compute CHC
1414 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1415 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 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 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1422 CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
1423 CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
1424 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 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 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1430 CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
1431 CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
1432 CALL cp_fm_to_fm_submat(matrix_mmt_fm, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
1433
1434 CALL cp_fm_release(matrix_mmt_fm)
1435
1436 ! (reuse matrix_sc and matrix_hc to computr HZ and SZ)
1437 CALL dbcsr_get_info(matrix_pz, nfullrows_total=n, nfullcols_total=k)
1438 CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, matrix_pz, 0.0_dp, matrix_hc, last_column=k)
1439 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 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 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1445 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 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 CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
1451 CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
1452
1453 CALL dbcsr_release_p(matrix_mm)
1454 CALL dbcsr_release_p(matrix_sc)
1455 CALL dbcsr_release_p(matrix_hc)
1456
1457 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 para_env=mo_coeff%matrix_struct%para_env)
1463 CALL cp_fm_create(matrix_nm_fm, fm_struct_tmp, name="m_nxm")
1464 CALL cp_fm_create(matrix_z_fm, fm_struct_tmp, name="m_nxm")
1465 CALL cp_fm_struct_release(fm_struct_tmp)
1466
1467 CALL copy_dbcsr_to_fm(matrix_pz, matrix_z_fm)
1468 ! extract egenvectors
1469 CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1, 1, 1, 1)
1470 CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, mo_notconv_fm, matrix_mm_fm, 0.0_dp, matrix_nm_fm)
1471 CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1 + nmat, 1, 1, 1)
1472 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 CALL dbcsr_release_p(matrix_z)
1475 CALL dbcsr_release_p(matrix_pz)
1476 CALL cp_fm_release(matrix_z_fm)
1477 CALL cp_fm_release(s_block)
1478 CALL cp_fm_release(h_block)
1479 CALL cp_fm_release(w_block)
1480 CALL cp_fm_release(v_block)
1481 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 IF (nmo_converged > 0) THEN
1485 jj = 1
1486 DO j = 1, nset_not_conv
1487 i_first = inotconv_set(j, 1)
1488 i_last = inotconv_set(j, 2)
1489 n = i_last - i_first + 1
1490 CALL cp_fm_to_fm_submat(matrix_nm_fm, mo_coeff, nao, n, 1, jj, 1, i_first)
1491 eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
1492 jj = jj + n
1493 END DO
1494 DEALLOCATE (iconv_set)
1495 DEALLOCATE (inotconv_set)
1496
1497 CALL dbcsr_release_p(mo_notconv)
1498 CALL dbcsr_release_p(c_out)
1499 CALL cp_fm_release(mo_notconv_fm)
1500 DEALLOCATE (mo_notconv_fm)
1501 ELSE
1502 CALL cp_fm_to_fm(matrix_nm_fm, mo_coeff)
1503 eigenvalues(1:nmo) = evals(1:nmo)
1504 END IF
1505 DEALLOCATE (evals)
1506
1507 CALL cp_fm_release(matrix_nm_fm)
1508 CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
1509
1510 t2 = m_walltime()
1511 IF (output_unit > 0) THEN
1512 WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1513 iteration, nmo_converged, max_norm, min_norm, t2 - t1
1514 END IF
1515 t1 = m_walltime()
1516
1517 END DO ! iteration
1518
1519 DEALLOCATE (ritz_coeff)
1520 DEALLOCATE (iconv)
1521 DEALLOCATE (inotconv)
1522 DEALLOCATE (vnorm)
1523
1524 CALL timestop(handle)
1525
1526 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 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 CALL timeset(routinen, handle)
1550
1551 CALL cp_fm_to_fm(s_block, w_block)
1552 CALL cp_fm_cholesky_decompose(s_block, info_out=info)
1553 IF (info == 0) THEN
1554 CALL cp_fm_triangular_invert(s_block)
1555 CALL cp_fm_cholesky_restore(h_block, ndim, s_block, w_block, "MULTIPLY", pos="RIGHT")
1556 CALL cp_fm_cholesky_restore(w_block, ndim, s_block, h_block, "MULTIPLY", pos="LEFT", transa="T")
1557 CALL choose_eigv_solver(h_block, w_block, evals)
1558 CALL cp_fm_cholesky_restore(w_block, ndim, s_block, v_block, "MULTIPLY")
1559 ELSE
1560! S^(-1/2)
1561 CALL cp_fm_power(w_block, s_block, -0.5_dp, 1.0e-5_dp, info)
1562 CALL cp_fm_to_fm(w_block, s_block)
1563 CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, h_block, s_block, 0.0_dp, w_block)
1564 CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, h_block)
1565 CALL choose_eigv_solver(h_block, w_block, evals)
1566 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 CALL timestop(handle)
1570
1571 END SUBROUTINE reduce_extended_space
1572
1573END MODULE qs_scf_block_davidson
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
preconditioner::apply_preconditioner
Definition preconditioner.F:115
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:243
cp_cfm_basic_linalg::cp_cfm_get_diag
subroutine, public cp_cfm_get_diag(matrix, diag)
returns the diagonal of a complex full matrix: diag(i)= A_{ii}. Each diagonal entry is owned by one p...
Definition cp_cfm_basic_linalg.F:1479
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:562
cp_cfm_basic_linalg::cp_cfm_vectorsnorm
subroutine, public cp_cfm_vectorsnorm(matrix, norm_array)
find the norm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*conjg(A_{ij}) ) Complex-valued mirror of...
Definition cp_cfm_basic_linalg.F:1412
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:632
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:340
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:357
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_set_submatrix
subroutine, public cp_cfm_set_submatrix(matrix, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
Set a sub-matrix of the full matrix: matrix(start_row:start_row+n_rows,start_col:start_col+n_cols) = ...
Definition cp_cfm_types.F:495
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:203
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_get_diag
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
Definition cp_dbcsr_contrib.F:544
cp_dbcsr_contrib::dbcsr_scale_by_vector
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
Definition cp_dbcsr_contrib.F:486
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:799
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:232
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_row_template
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
Definition cp_dbcsr_operations.F:1186
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:1143
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1532
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition cp_fm_basic_linalg.F:1968
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition cp_fm_basic_linalg.F:568
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:313
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1144
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_vectorsnorm
subroutine, public cp_fm_vectorsnorm(matrix, norm_array)
find the inorm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
Definition cp_fm_types.F:1213
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1455
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
qs_block_davidson_types
module that contains the algorithms to perform an iterative diagonalization by the block-Davidson app...
Definition qs_block_davidson_types.F:18
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
qs_scf_block_davidson
module that contains the algorithms to perform an iterative diagonalization by the block-Davidson app...
Definition qs_scf_block_davidson.F:18
qs_scf_block_davidson::generate_extended_space_c
subroutine, public generate_extended_space_c(bdav_env, mos, matrix_h, matrix_s, output_unit, eps_iter, eps_iter_empty, preconditioner)
iterative diagonalization by the block-Davidson approach for one complex K point; complex counterpart...
Definition qs_scf_block_davidson.F:501
qs_scf_block_davidson::generate_extended_space_sparse
subroutine, public generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:1069
qs_scf_block_davidson::generate_extended_space
subroutine, public generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition qs_scf_block_davidson.F:98
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
message_passing::mp_comm_type
Definition message_passing.F:179
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:60
qs_block_davidson_types::davidson_type
Definition qs_block_davidson_types.F:38
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53