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_dbcsr_api, ONLY: &
21 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_init_p, &
22 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
23 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release_p, dbcsr_type, &
24 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
25 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_diag,&
26 : dbcsr_scale_by_vector
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr,&
29 : cp_dbcsr_m_by_n_from_row_template,&
30 : cp_dbcsr_m_by_n_from_template,&
31 : cp_dbcsr_sm_fm_multiply
32 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
33 : cp_fm_scale_and_add,&
34 : cp_fm_symm,&
35 : cp_fm_transpose,&
36 : cp_fm_triangular_invert,&
37 : cp_fm_uplo_to_full
38 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
39 : cp_fm_cholesky_restore
40 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
41 : cp_fm_power
42 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 : cp_fm_struct_release,&
44 : cp_fm_struct_type
45 : USE cp_fm_types, ONLY: cp_fm_create,&
46 : cp_fm_get_diag,&
47 : cp_fm_release,&
48 : cp_fm_set_all,&
49 : cp_fm_to_fm,&
50 : cp_fm_to_fm_submat,&
51 : cp_fm_type,&
52 : cp_fm_vectorsnorm
53 : USE kinds, ONLY: dp
54 : USE machine, ONLY: m_walltime
55 : USE message_passing, ONLY: mp_comm_type
56 : USE parallel_gemm_api, ONLY: parallel_gemm
57 : USE preconditioner, ONLY: apply_preconditioner
58 : USE preconditioner_types, ONLY: preconditioner_type
59 : USE qs_block_davidson_types, ONLY: davidson_type
60 : USE qs_mo_types, ONLY: get_mo_set,&
61 : mo_set_type
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 : PRIVATE
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_block_davidson'
67 :
68 : PUBLIC :: generate_extended_space, generate_extended_space_sparse
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param bdav_env ...
75 : !> \param mo_set ...
76 : !> \param matrix_h ...
77 : !> \param matrix_s ...
78 : !> \param output_unit ...
79 : !> \param preconditioner ...
80 : ! **************************************************************************************************
81 40 : SUBROUTINE generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
82 : preconditioner)
83 :
84 : TYPE(davidson_type) :: bdav_env
85 : TYPE(mo_set_type), INTENT(IN) :: mo_set
86 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
87 : INTEGER, INTENT(IN) :: output_unit
88 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
89 :
90 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space'
91 :
92 : INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, max_iter, n, nao, nmat, nmat2, &
93 : nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
94 40 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
95 40 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
96 : LOGICAL :: converged, do_apply_preconditioner
97 : REAL(dp) :: lambda, max_norm, min_norm, t1, t2
98 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm
99 40 : REAL(dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals
100 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
101 : TYPE(cp_fm_type) :: c_conv, c_notconv, c_out, h_block, h_fm, &
102 : m_hc, m_sc, m_tmp, mt_tmp, s_block, &
103 : s_fm, v_block, w_block
104 : TYPE(cp_fm_type), POINTER :: c_pz, c_z, mo_coeff
105 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
106 :
107 40 : CALL timeset(routineN, handle)
108 :
109 40 : NULLIFY (mo_coeff, mo_coeff_b, eigenvalues)
110 :
111 40 : do_apply_preconditioner = .FALSE.
112 40 : IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE.
113 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, eigenvalues=eigenvalues, &
114 40 : nao=nao, nmo=nmo, homo=homo)
115 40 : IF (do_apply_preconditioner) THEN
116 36 : max_iter = bdav_env%max_iter
117 : ELSE
118 : max_iter = 1
119 : END IF
120 :
121 40 : NULLIFY (c_z, c_pz)
122 40 : NULLIFY (evals, eig_not_conv)
123 40 : t1 = m_walltime()
124 40 : IF (output_unit > 0) THEN
125 : WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
126 0 : " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
127 : END IF
128 :
129 120 : ALLOCATE (iconv(nmo))
130 80 : ALLOCATE (inotconv(nmo))
131 120 : ALLOCATE (ritz_coeff(nmo))
132 80 : ALLOCATE (vnorm(nmo))
133 :
134 40 : converged = .FALSE.
135 124 : DO iter = 1, max_iter
136 :
137 : ! compute Ritz values
138 88 : ritz_coeff = 0.0_dp
139 88 : CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc")
140 88 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo)
141 88 : CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc")
142 88 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_coeff, m_sc, nmo)
143 :
144 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
145 : context=mo_coeff%matrix_struct%context, &
146 88 : para_env=mo_coeff%matrix_struct%para_env)
147 88 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
148 88 : CALL cp_fm_struct_release(fm_struct_tmp)
149 :
150 88 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp)
151 88 : CALL cp_fm_get_diag(m_tmp, ritz_coeff)
152 88 : CALL cp_fm_release(m_tmp)
153 :
154 : ! Check for converged eigenvectors
155 88 : c_z => bdav_env%matrix_z
156 88 : c_pz => bdav_env%matrix_pz
157 88 : CALL cp_fm_to_fm(m_sc, c_z)
158 88 : CALL cp_fm_column_scale(c_z, ritz_coeff)
159 88 : CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc)
160 88 : CALL cp_fm_vectorsnorm(c_z, vnorm)
161 :
162 88 : nmo_converged = 0
163 88 : nmo_not_converged = 0
164 88 : max_norm = 0.0_dp
165 88 : min_norm = 1.e10_dp
166 3256 : DO imo = 1, nmo
167 3168 : max_norm = MAX(max_norm, vnorm(imo))
168 3256 : min_norm = MIN(min_norm, vnorm(imo))
169 : END DO
170 88 : iconv = 0
171 88 : inotconv = 0
172 3256 : DO imo = 1, nmo
173 3256 : IF (vnorm(imo) <= bdav_env%eps_iter) THEN
174 136 : nmo_converged = nmo_converged + 1
175 136 : iconv(nmo_converged) = imo
176 : ELSE
177 3032 : nmo_not_converged = nmo_not_converged + 1
178 3032 : inotconv(nmo_not_converged) = imo
179 : END IF
180 : END DO
181 :
182 88 : IF (nmo_converged > 0) THEN
183 48 : ALLOCATE (iconv_set(nmo_converged, 2))
184 48 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
185 16 : i_last = iconv(1)
186 16 : nset = 0
187 152 : DO j = 1, nmo_converged
188 136 : imo = iconv(j)
189 :
190 152 : IF (imo == i_last + 1) THEN
191 102 : i_last = imo
192 102 : iconv_set(nset, 2) = imo
193 : ELSE
194 34 : i_last = imo
195 34 : nset = nset + 1
196 34 : iconv_set(nset, 1) = imo
197 34 : iconv_set(nset, 2) = imo
198 : END IF
199 : END DO
200 16 : nset_conv = nset
201 :
202 16 : i_last = inotconv(1)
203 16 : nset = 0
204 456 : DO j = 1, nmo_not_converged
205 440 : imo = inotconv(j)
206 :
207 456 : IF (imo == i_last + 1) THEN
208 398 : i_last = imo
209 398 : inotconv_set(nset, 2) = imo
210 : ELSE
211 42 : i_last = imo
212 42 : nset = nset + 1
213 42 : inotconv_set(nset, 1) = imo
214 42 : inotconv_set(nset, 2) = imo
215 : END IF
216 : END DO
217 16 : nset_not_conv = nset
218 16 : CALL cp_fm_release(m_sc)
219 16 : CALL cp_fm_release(m_hc)
220 16 : NULLIFY (c_z, c_pz)
221 : END IF
222 :
223 88 : IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
224 4 : converged = .TRUE.
225 4 : DEALLOCATE (iconv_set)
226 4 : DEALLOCATE (inotconv_set)
227 4 : t2 = m_walltime()
228 4 : IF (output_unit > 0) THEN
229 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
230 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
231 :
232 0 : WRITE (output_unit, *) " Reached convergence in ", iter, &
233 0 : " Davidson iterations"
234 : END IF
235 :
236 : EXIT
237 : END IF
238 :
239 84 : IF (nmo_converged > 0) THEN
240 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
241 : context=mo_coeff%matrix_struct%context, &
242 12 : para_env=mo_coeff%matrix_struct%para_env)
243 : !allocate h_fm
244 12 : CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp")
245 : !allocate s_fm
246 12 : CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp")
247 : !copy matrix_h in h_fm
248 12 : CALL copy_dbcsr_to_fm(matrix_h, h_fm)
249 12 : CALL cp_fm_uplo_to_full(h_fm, s_fm)
250 :
251 : !copy matrix_s in s_fm
252 : ! CALL cp_fm_set_all(s_fm,0.0_dp)
253 12 : CALL copy_dbcsr_to_fm(matrix_s, s_fm)
254 :
255 : !allocate c_out
256 12 : CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp")
257 : ! set c_out to zero
258 12 : CALL cp_fm_set_all(c_out, 0.0_dp)
259 12 : CALL cp_fm_struct_release(fm_struct_tmp)
260 :
261 : !allocate c_conv
262 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
263 : context=mo_coeff%matrix_struct%context, &
264 12 : para_env=mo_coeff%matrix_struct%para_env)
265 12 : CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv")
266 12 : CALL cp_fm_set_all(c_conv, 0.0_dp)
267 : !allocate m_tmp
268 12 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc")
269 12 : CALL cp_fm_struct_release(fm_struct_tmp)
270 : END IF
271 :
272 : !allocate c_notconv
273 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
274 : context=mo_coeff%matrix_struct%context, &
275 84 : para_env=mo_coeff%matrix_struct%para_env)
276 84 : CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv")
277 84 : CALL cp_fm_set_all(c_notconv, 0.0_dp)
278 84 : IF (nmo_converged > 0) THEN
279 12 : CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc")
280 12 : CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc")
281 : !allocate c_z
282 12 : ALLOCATE (c_z, c_pz)
283 12 : CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z")
284 12 : CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz")
285 12 : CALL cp_fm_set_all(c_z, 0.0_dp)
286 :
287 : ! sum contributions to c_out
288 12 : jj = 1
289 34 : DO j = 1, nset_conv
290 22 : i_first = iconv_set(j, 1)
291 22 : i_last = iconv_set(j, 2)
292 22 : n = i_last - i_first + 1
293 22 : CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj)
294 34 : jj = jj + n
295 : END DO
296 12 : CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp)
297 12 : CALL parallel_gemm('N', 'T', nao, nao, nmo_converged, 1.0_dp, m_tmp, m_tmp, 0.0_dp, c_out)
298 :
299 : ! project c_out out of H
300 12 : lambda = 100.0_dp*ABS(eigenvalues(homo))
301 12 : CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm)
302 12 : CALL cp_fm_release(m_tmp)
303 12 : CALL cp_fm_release(h_fm)
304 :
305 : END IF
306 :
307 : !allocate m_tmp
308 84 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm")
309 84 : CALL cp_fm_struct_release(fm_struct_tmp)
310 84 : IF (nmo_converged > 0) THEN
311 36 : ALLOCATE (eig_not_conv(nmo_not_converged))
312 12 : jj = 1
313 42 : DO j = 1, nset_not_conv
314 30 : i_first = inotconv_set(j, 1)
315 30 : i_last = inotconv_set(j, 2)
316 30 : n = i_last - i_first + 1
317 30 : CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj)
318 404 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
319 42 : jj = jj + n
320 : END DO
321 12 : CALL parallel_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc)
322 12 : CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_notconv, 0.0_dp, m_sc)
323 : ! extend suspace using only the not converged vectors
324 12 : CALL cp_fm_to_fm(m_sc, m_tmp)
325 12 : CALL cp_fm_column_scale(m_tmp, eig_not_conv)
326 12 : CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc)
327 12 : DEALLOCATE (eig_not_conv)
328 12 : CALL cp_fm_to_fm(m_tmp, c_z)
329 : ELSE
330 72 : CALL cp_fm_to_fm(mo_coeff, c_notconv)
331 : END IF
332 :
333 : !preconditioner
334 84 : IF (do_apply_preconditioner) THEN
335 80 : IF (preconditioner%in_use /= 0) THEN
336 80 : CALL apply_preconditioner(preconditioner, c_z, c_pz)
337 : ELSE
338 0 : CALL cp_fm_to_fm(c_z, c_pz)
339 : END IF
340 : ELSE
341 4 : CALL cp_fm_to_fm(c_z, c_pz)
342 : END IF
343 84 : CALL cp_fm_release(m_tmp)
344 :
345 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo_not_converged, ncol_global=nmo_not_converged, &
346 : context=mo_coeff%matrix_struct%context, &
347 84 : para_env=mo_coeff%matrix_struct%para_env)
348 :
349 84 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm")
350 84 : CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm")
351 84 : CALL cp_fm_struct_release(fm_struct_tmp)
352 :
353 84 : nmat = nmo_not_converged
354 84 : nmat2 = 2*nmo_not_converged
355 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
356 : context=mo_coeff%matrix_struct%context, &
357 84 : para_env=mo_coeff%matrix_struct%para_env)
358 :
359 84 : CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
360 84 : CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
361 84 : CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
362 84 : CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
363 252 : ALLOCATE (evals(nmat2))
364 :
365 84 : CALL cp_fm_struct_release(fm_struct_tmp)
366 :
367 : ! compute CSC
368 84 : CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
369 :
370 : ! compute CHC
371 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp)
372 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1)
373 :
374 : ! compute ZSC
375 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
376 84 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
377 84 : CALL cp_fm_transpose(m_tmp, mt_tmp)
378 84 : CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
379 : ! compute ZHC
380 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
381 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
382 84 : CALL cp_fm_transpose(m_tmp, mt_tmp)
383 84 : CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
384 :
385 84 : CALL cp_fm_release(mt_tmp)
386 :
387 : ! reuse m_sc and m_hc to computr HZ and SZ
388 84 : IF (nmo_converged > 0) THEN
389 12 : CALL parallel_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc)
390 12 : CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc)
391 :
392 12 : CALL cp_fm_release(c_out)
393 12 : CALL cp_fm_release(c_conv)
394 12 : CALL cp_fm_release(s_fm)
395 : ELSE
396 72 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo)
397 72 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo)
398 : END IF
399 :
400 : ! compute ZSZ
401 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
402 84 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
403 : ! compute ZHZ
404 84 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
405 84 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
406 :
407 84 : CALL cp_fm_release(m_sc)
408 :
409 : ! solution of the reduced eigenvalues problem
410 84 : CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
411 :
412 : ! extract egenvectors
413 84 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1)
414 84 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc)
415 84 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1)
416 84 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc)
417 :
418 84 : CALL cp_fm_release(m_tmp)
419 :
420 84 : CALL cp_fm_release(c_notconv)
421 84 : CALL cp_fm_release(s_block)
422 84 : CALL cp_fm_release(h_block)
423 84 : CALL cp_fm_release(w_block)
424 84 : CALL cp_fm_release(v_block)
425 :
426 84 : IF (nmo_converged > 0) THEN
427 12 : CALL cp_fm_release(c_z)
428 12 : CALL cp_fm_release(c_pz)
429 12 : DEALLOCATE (c_z, c_pz)
430 12 : jj = 1
431 42 : DO j = 1, nset_not_conv
432 30 : i_first = inotconv_set(j, 1)
433 30 : i_last = inotconv_set(j, 2)
434 30 : n = i_last - i_first + 1
435 30 : CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first)
436 808 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
437 42 : jj = jj + n
438 : END DO
439 12 : DEALLOCATE (iconv_set)
440 12 : DEALLOCATE (inotconv_set)
441 : ELSE
442 72 : CALL cp_fm_to_fm(m_hc, mo_coeff)
443 5328 : eigenvalues(1:nmo) = evals(1:nmo)
444 : END IF
445 84 : DEALLOCATE (evals)
446 :
447 84 : CALL cp_fm_release(m_hc)
448 :
449 84 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
450 :
451 84 : t2 = m_walltime()
452 84 : IF (output_unit > 0) THEN
453 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
454 0 : iter, nmo_converged, max_norm, min_norm, t2 - t1
455 : END IF
456 720 : t1 = m_walltime()
457 :
458 : END DO ! iter
459 :
460 40 : DEALLOCATE (iconv)
461 40 : DEALLOCATE (inotconv)
462 40 : DEALLOCATE (ritz_coeff)
463 40 : DEALLOCATE (vnorm)
464 :
465 40 : CALL timestop(handle)
466 120 : END SUBROUTINE generate_extended_space
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param bdav_env ...
471 : !> \param mo_set ...
472 : !> \param matrix_h ...
473 : !> \param matrix_s ...
474 : !> \param output_unit ...
475 : !> \param preconditioner ...
476 : ! **************************************************************************************************
477 64 : SUBROUTINE generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
478 : preconditioner)
479 :
480 : TYPE(davidson_type) :: bdav_env
481 : TYPE(mo_set_type), INTENT(IN) :: mo_set
482 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
483 : INTEGER, INTENT(IN) :: output_unit
484 : TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
485 :
486 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_extended_space_sparse'
487 :
488 : INTEGER :: col_offset, handle, homo, i_first, i_last, imo, iteration, j, jj, k, max_iter, n, &
489 : nao, nmat, nmat2, nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
490 64 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
491 64 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
492 : LOGICAL :: converged, do_apply_preconditioner
493 : REAL(dp) :: lambda, max_norm, min_norm, t1, t2
494 64 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, ritz_coeff, vnorm
495 64 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues
496 64 : REAL(dp), DIMENSION(:, :), POINTER :: block
497 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
498 : TYPE(cp_fm_type) :: h_block, matrix_mm_fm, matrix_mmt_fm, &
499 : matrix_nm_fm, matrix_z_fm, mo_conv_fm, &
500 : s_block, v_block, w_block
501 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_notconv_fm
502 : TYPE(dbcsr_iterator_type) :: iter
503 : TYPE(dbcsr_type), POINTER :: c_out, matrix_hc, matrix_mm, matrix_pz, &
504 : matrix_sc, matrix_z, mo_coeff_b, &
505 : mo_conv, mo_notconv, smo_conv
506 : TYPE(mp_comm_type) :: group
507 :
508 64 : CALL timeset(routineN, handle)
509 :
510 64 : do_apply_preconditioner = .FALSE.
511 64 : IF (PRESENT(preconditioner)) do_apply_preconditioner = .TRUE.
512 :
513 64 : NULLIFY (mo_coeff, mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm)
514 64 : NULLIFY (mo_notconv_fm, mo_conv, mo_notconv, smo_conv, c_out)
515 64 : NULLIFY (fm_struct_tmp)
516 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
517 64 : eigenvalues=eigenvalues, homo=homo, nao=nao, nmo=nmo)
518 64 : IF (do_apply_preconditioner) THEN
519 56 : max_iter = bdav_env%max_iter
520 : ELSE
521 : max_iter = 1
522 : END IF
523 :
524 64 : t1 = m_walltime()
525 64 : IF (output_unit > 0) THEN
526 : WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
527 0 : " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", REPEAT("-", 60)
528 : END IF
529 :
530 : ! Allocate array for Ritz values
531 192 : ALLOCATE (ritz_coeff(nmo))
532 192 : ALLOCATE (iconv(nmo))
533 128 : ALLOCATE (inotconv(nmo))
534 128 : ALLOCATE (vnorm(nmo))
535 :
536 64 : converged = .FALSE.
537 128 : DO iteration = 1, max_iter
538 88 : NULLIFY (c_out, mo_conv, mo_notconv_fm, mo_notconv)
539 : ! Prepare HC and SC, using mo_coeff_b (sparse), these are still sparse
540 88 : CALL dbcsr_init_p(matrix_hc)
541 : CALL dbcsr_create(matrix_hc, template=mo_coeff_b, &
542 : name="matrix_hc", &
543 88 : matrix_type=dbcsr_type_no_symmetry)
544 88 : CALL dbcsr_init_p(matrix_sc)
545 : CALL dbcsr_create(matrix_sc, template=mo_coeff_b, &
546 : name="matrix_sc", &
547 88 : matrix_type=dbcsr_type_no_symmetry)
548 :
549 88 : CALL dbcsr_get_info(mo_coeff_b, nfullrows_total=n, nfullcols_total=k, group=group)
550 88 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h, mo_coeff_b, 0.0_dp, matrix_hc, last_column=k)
551 88 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_coeff_b, 0.0_dp, matrix_sc, last_column=k)
552 :
553 : ! compute Ritz values
554 88 : ritz_coeff = 0.0_dp
555 : ! Allocate Sparse matrices: nmoxnmo
556 : ! matrix_mm
557 :
558 88 : CALL dbcsr_init_p(matrix_mm)
559 : CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo, n=nmo, &
560 88 : sym=dbcsr_type_no_symmetry)
561 :
562 88 : CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_coeff_b, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
563 88 : CALL dbcsr_get_diag(matrix_mm, ritz_coeff)
564 88 : CALL mo_coeff%matrix_struct%para_env%sum(ritz_coeff)
565 :
566 : ! extended subspace P Z = P [H - theta S]C this ia another matrix of type and size as mo_coeff_b
567 88 : CALL dbcsr_init_p(matrix_z)
568 : CALL dbcsr_create(matrix_z, template=mo_coeff_b, &
569 : name="matrix_z", &
570 88 : matrix_type=dbcsr_type_no_symmetry)
571 88 : CALL dbcsr_copy(matrix_z, matrix_sc)
572 88 : CALL dbcsr_scale_by_vector(matrix_z, ritz_coeff, side='right')
573 88 : CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
574 :
575 : ! Compute the column norms of matrix_z.
576 88 : vnorm = 0.0_dp
577 88 : CALL dbcsr_iterator_start(iter, matrix_z)
578 792 : DO WHILE (dbcsr_iterator_blocks_left(iter))
579 704 : CALL dbcsr_iterator_next_block(iter, block=block, col_offset=col_offset)
580 13464 : DO j = 1, SIZE(block, 2)
581 178112 : vnorm(col_offset + j - 1) = vnorm(col_offset + j - 1) + SUM(block(:, j)**2)
582 : END DO
583 : END DO
584 88 : CALL dbcsr_iterator_stop(iter)
585 88 : CALL group%sum(vnorm)
586 3256 : vnorm = SQRT(vnorm)
587 :
588 : ! Check for converged eigenvectors
589 88 : nmo_converged = 0
590 88 : nmo_not_converged = 0
591 88 : max_norm = 0.0_dp
592 88 : min_norm = 1.e10_dp
593 3256 : DO imo = 1, nmo
594 3168 : max_norm = MAX(max_norm, vnorm(imo))
595 3256 : min_norm = MIN(min_norm, vnorm(imo))
596 : END DO
597 88 : iconv = 0
598 88 : inotconv = 0
599 :
600 3256 : DO imo = 1, nmo
601 3256 : IF (vnorm(imo) <= bdav_env%eps_iter) THEN
602 836 : nmo_converged = nmo_converged + 1
603 836 : iconv(nmo_converged) = imo
604 : ELSE
605 2332 : nmo_not_converged = nmo_not_converged + 1
606 2332 : inotconv(nmo_not_converged) = imo
607 : END IF
608 : END DO
609 :
610 88 : IF (nmo_converged > 0) THEN
611 90 : ALLOCATE (iconv_set(nmo_converged, 2))
612 88 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
613 30 : i_last = iconv(1)
614 30 : nset = 0
615 866 : DO j = 1, nmo_converged
616 836 : imo = iconv(j)
617 :
618 866 : IF (imo == i_last + 1) THEN
619 772 : i_last = imo
620 772 : iconv_set(nset, 2) = imo
621 : ELSE
622 64 : i_last = imo
623 64 : nset = nset + 1
624 64 : iconv_set(nset, 1) = imo
625 64 : iconv_set(nset, 2) = imo
626 : END IF
627 : END DO
628 30 : nset_conv = nset
629 :
630 30 : i_last = inotconv(1)
631 30 : nset = 0
632 274 : DO j = 1, nmo_not_converged
633 244 : imo = inotconv(j)
634 :
635 274 : IF (imo == i_last + 1) THEN
636 184 : i_last = imo
637 184 : inotconv_set(nset, 2) = imo
638 : ELSE
639 60 : i_last = imo
640 60 : nset = nset + 1
641 60 : inotconv_set(nset, 1) = imo
642 60 : inotconv_set(nset, 2) = imo
643 : END IF
644 : END DO
645 30 : nset_not_conv = nset
646 :
647 30 : CALL dbcsr_release_p(matrix_hc)
648 30 : CALL dbcsr_release_p(matrix_sc)
649 30 : CALL dbcsr_release_p(matrix_z)
650 30 : CALL dbcsr_release_p(matrix_mm)
651 : END IF
652 :
653 88 : IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
654 24 : DEALLOCATE (iconv_set)
655 :
656 24 : DEALLOCATE (inotconv_set)
657 :
658 24 : converged = .TRUE.
659 24 : t2 = m_walltime()
660 24 : IF (output_unit > 0) THEN
661 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
662 0 : iteration, nmo_converged, max_norm, min_norm, t2 - t1
663 :
664 0 : WRITE (output_unit, *) " Reached convergence in ", iteration, &
665 0 : " Davidson iterations"
666 : END IF
667 :
668 : EXIT
669 : END IF
670 :
671 64 : IF (nmo_converged > 0) THEN
672 :
673 : !allocate mo_conv_fm
674 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
675 : context=mo_coeff%matrix_struct%context, &
676 6 : para_env=mo_coeff%matrix_struct%para_env)
677 6 : CALL cp_fm_create(mo_conv_fm, fm_struct_tmp, name="mo_conv_fm")
678 :
679 6 : CALL cp_fm_struct_release(fm_struct_tmp)
680 :
681 : ! extract mo_conv from mo_coeff full matrix
682 6 : jj = 1
683 22 : DO j = 1, nset_conv
684 16 : i_first = iconv_set(j, 1)
685 16 : i_last = iconv_set(j, 2)
686 16 : n = i_last - i_first + 1
687 16 : CALL cp_fm_to_fm_submat(mo_coeff, mo_conv_fm, nao, n, 1, i_first, 1, jj)
688 22 : jj = jj + n
689 : END DO
690 :
691 : ! allocate c_out sparse matrix, to project out the converged MOS
692 6 : CALL dbcsr_init_p(c_out)
693 : CALL dbcsr_create(c_out, template=matrix_s, &
694 : name="c_out", &
695 6 : matrix_type=dbcsr_type_symmetric)
696 :
697 : ! allocate mo_conv sparse
698 6 : CALL dbcsr_init_p(mo_conv)
699 : CALL cp_dbcsr_m_by_n_from_row_template(mo_conv, template=matrix_s, n=nmo_converged, &
700 6 : sym=dbcsr_type_no_symmetry)
701 :
702 6 : CALL dbcsr_init_p(smo_conv)
703 : CALL cp_dbcsr_m_by_n_from_row_template(smo_conv, template=matrix_s, n=nmo_converged, &
704 6 : sym=dbcsr_type_no_symmetry)
705 :
706 6 : CALL copy_fm_to_dbcsr(mo_conv_fm, mo_conv) !fm->dbcsr
707 :
708 6 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_conv, 0.0_dp, smo_conv, last_column=nmo_converged)
709 6 : CALL dbcsr_multiply('n', 't', 1.0_dp, smo_conv, smo_conv, 0.0_dp, c_out, last_column=nao)
710 : ! project c_out out of H
711 6 : lambda = 100.0_dp*ABS(eigenvalues(homo))
712 6 : CALL dbcsr_add(c_out, matrix_h, lambda, 1.0_dp)
713 :
714 6 : CALL dbcsr_release_p(mo_conv)
715 6 : CALL dbcsr_release_p(smo_conv)
716 6 : CALL cp_fm_release(mo_conv_fm)
717 :
718 : !allocate c_notconv_fm
719 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
720 : context=mo_coeff%matrix_struct%context, &
721 6 : para_env=mo_coeff%matrix_struct%para_env)
722 6 : ALLOCATE (mo_notconv_fm)
723 6 : CALL cp_fm_create(mo_notconv_fm, fm_struct_tmp, name="mo_notconv_fm")
724 6 : CALL cp_fm_struct_release(fm_struct_tmp)
725 :
726 : ! extract mo_notconv from mo_coeff full matrix
727 18 : ALLOCATE (eig_not_conv(nmo_not_converged))
728 :
729 6 : jj = 1
730 24 : DO j = 1, nset_not_conv
731 18 : i_first = inotconv_set(j, 1)
732 18 : i_last = inotconv_set(j, 2)
733 18 : n = i_last - i_first + 1
734 18 : CALL cp_fm_to_fm_submat(mo_coeff, mo_notconv_fm, nao, n, 1, i_first, 1, jj)
735 186 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
736 24 : jj = jj + n
737 : END DO
738 :
739 : ! allocate mo_conv sparse
740 6 : CALL dbcsr_init_p(mo_notconv)
741 : CALL cp_dbcsr_m_by_n_from_row_template(mo_notconv, template=matrix_s, n=nmo_not_converged, &
742 6 : sym=dbcsr_type_no_symmetry)
743 :
744 6 : CALL dbcsr_init_p(matrix_hc)
745 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_hc, template=matrix_s, n=nmo_not_converged, &
746 6 : sym=dbcsr_type_no_symmetry)
747 :
748 6 : CALL dbcsr_init_p(matrix_sc)
749 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_sc, template=matrix_s, n=nmo_not_converged, &
750 6 : sym=dbcsr_type_no_symmetry)
751 :
752 6 : CALL dbcsr_init_p(matrix_z)
753 : CALL cp_dbcsr_m_by_n_from_row_template(matrix_z, template=matrix_s, n=nmo_not_converged, &
754 6 : sym=dbcsr_type_no_symmetry)
755 :
756 6 : CALL copy_fm_to_dbcsr(mo_notconv_fm, mo_notconv) !fm->dbcsr
757 :
758 : CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, mo_notconv, 0.0_dp, matrix_hc, &
759 6 : last_column=nmo_not_converged)
760 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_notconv, 0.0_dp, matrix_sc, &
761 6 : last_column=nmo_not_converged)
762 :
763 6 : CALL dbcsr_copy(matrix_z, matrix_sc)
764 6 : CALL dbcsr_scale_by_vector(matrix_z, eig_not_conv, side='right')
765 6 : CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
766 :
767 6 : DEALLOCATE (eig_not_conv)
768 :
769 : ! matrix_mm
770 6 : CALL dbcsr_init_p(matrix_mm)
771 : CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo_not_converged, n=nmo_not_converged, &
772 6 : sym=dbcsr_type_no_symmetry)
773 :
774 : CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_notconv, matrix_hc, 0.0_dp, matrix_mm, &
775 18 : last_column=nmo_not_converged)
776 :
777 : ELSE
778 58 : mo_notconv => mo_coeff_b
779 58 : mo_notconv_fm => mo_coeff
780 58 : c_out => matrix_h
781 : END IF
782 :
783 : ! allocate matrix_pz using as template matrix_z
784 64 : CALL dbcsr_init_p(matrix_pz)
785 : CALL dbcsr_create(matrix_pz, template=matrix_z, &
786 : name="matrix_pz", &
787 64 : matrix_type=dbcsr_type_no_symmetry)
788 :
789 64 : IF (do_apply_preconditioner) THEN
790 56 : IF (preconditioner%in_use /= 0) THEN
791 56 : CALL apply_preconditioner(preconditioner, matrix_z, matrix_pz)
792 : ELSE
793 0 : CALL dbcsr_copy(matrix_pz, matrix_z)
794 : END IF
795 : ELSE
796 8 : CALL dbcsr_copy(matrix_pz, matrix_z)
797 : END IF
798 :
799 : !allocate NMOxNMO full matrices
800 64 : nmat = nmo_not_converged
801 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat, ncol_global=nmat, &
802 : context=mo_coeff%matrix_struct%context, &
803 64 : para_env=mo_coeff%matrix_struct%para_env)
804 64 : CALL cp_fm_create(matrix_mm_fm, fm_struct_tmp, name="m_tmp_mxm")
805 64 : CALL cp_fm_create(matrix_mmt_fm, fm_struct_tmp, name="mt_tmp_mxm")
806 64 : CALL cp_fm_struct_release(fm_struct_tmp)
807 :
808 : !allocate 2NMOx2NMO full matrices
809 64 : nmat2 = 2*nmo_not_converged
810 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
811 : context=mo_coeff%matrix_struct%context, &
812 64 : para_env=mo_coeff%matrix_struct%para_env)
813 :
814 64 : CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
815 64 : CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
816 64 : CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
817 64 : CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
818 192 : ALLOCATE (evals(nmat2))
819 64 : CALL cp_fm_struct_release(fm_struct_tmp)
820 :
821 : ! compute CSC
822 64 : CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
823 : ! compute CHC
824 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
825 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1, 1)
826 :
827 : ! compute the bottom left ZSC (top right is transpose)
828 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=nmat)
829 : ! set the bottom left part of S[C,Z] block matrix ZSC
830 : !copy sparse to full
831 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
832 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
833 64 : CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
834 64 : CALL cp_fm_to_fm_submat(matrix_mmt_fm, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
835 :
836 : ! compute the bottom left ZHC (top right is transpose)
837 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=nmat)
838 : ! set the bottom left part of S[C,Z] block matrix ZHC
839 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
840 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
841 64 : CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
842 64 : CALL cp_fm_to_fm_submat(matrix_mmt_fm, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
843 :
844 64 : CALL cp_fm_release(matrix_mmt_fm)
845 :
846 : ! (reuse matrix_sc and matrix_hc to computr HZ and SZ)
847 64 : CALL dbcsr_get_info(matrix_pz, nfullrows_total=n, nfullcols_total=k)
848 64 : CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, matrix_pz, 0.0_dp, matrix_hc, last_column=k)
849 64 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, matrix_pz, 0.0_dp, matrix_sc, last_column=k)
850 :
851 : ! compute the bottom right ZSZ
852 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=k)
853 : ! set the bottom right part of S[C,Z] block matrix ZSZ
854 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
855 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
856 :
857 : ! compute the bottom right ZHZ
858 64 : CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
859 : ! set the bottom right part of H[C,Z] block matrix ZHZ
860 64 : CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
861 64 : CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
862 :
863 64 : CALL dbcsr_release_p(matrix_mm)
864 64 : CALL dbcsr_release_p(matrix_sc)
865 64 : CALL dbcsr_release_p(matrix_hc)
866 :
867 64 : CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
868 :
869 : ! allocate two (nao x nmat) full matrix
870 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmat, &
871 : context=mo_coeff%matrix_struct%context, &
872 64 : para_env=mo_coeff%matrix_struct%para_env)
873 64 : CALL cp_fm_create(matrix_nm_fm, fm_struct_tmp, name="m_nxm")
874 64 : CALL cp_fm_create(matrix_z_fm, fm_struct_tmp, name="m_nxm")
875 64 : CALL cp_fm_struct_release(fm_struct_tmp)
876 :
877 64 : CALL copy_dbcsr_to_fm(matrix_pz, matrix_z_fm)
878 : ! extract egenvectors
879 64 : CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1, 1, 1, 1)
880 64 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, mo_notconv_fm, matrix_mm_fm, 0.0_dp, matrix_nm_fm)
881 64 : CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1 + nmat, 1, 1, 1)
882 64 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, matrix_z_fm, matrix_mm_fm, 1.0_dp, matrix_nm_fm)
883 :
884 64 : CALL dbcsr_release_p(matrix_z)
885 64 : CALL dbcsr_release_p(matrix_pz)
886 64 : CALL cp_fm_release(matrix_z_fm)
887 64 : CALL cp_fm_release(s_block)
888 64 : CALL cp_fm_release(h_block)
889 64 : CALL cp_fm_release(w_block)
890 64 : CALL cp_fm_release(v_block)
891 64 : CALL cp_fm_release(matrix_mm_fm)
892 :
893 : ! in case some vector are already converged only a subset of vectors are copied in the MOS
894 64 : IF (nmo_converged > 0) THEN
895 6 : jj = 1
896 24 : DO j = 1, nset_not_conv
897 18 : i_first = inotconv_set(j, 1)
898 18 : i_last = inotconv_set(j, 2)
899 18 : n = i_last - i_first + 1
900 18 : CALL cp_fm_to_fm_submat(matrix_nm_fm, mo_coeff, nao, n, 1, jj, 1, i_first)
901 186 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
902 24 : jj = jj + n
903 : END DO
904 6 : DEALLOCATE (iconv_set)
905 6 : DEALLOCATE (inotconv_set)
906 :
907 6 : CALL dbcsr_release_p(mo_notconv)
908 6 : CALL dbcsr_release_p(c_out)
909 6 : CALL cp_fm_release(mo_notconv_fm)
910 6 : DEALLOCATE (mo_notconv_fm)
911 : ELSE
912 58 : CALL cp_fm_to_fm(matrix_nm_fm, mo_coeff)
913 2146 : eigenvalues(1:nmo) = evals(1:nmo)
914 : END IF
915 64 : DEALLOCATE (evals)
916 :
917 64 : CALL cp_fm_release(matrix_nm_fm)
918 64 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
919 :
920 64 : t2 = m_walltime()
921 64 : IF (output_unit > 0) THEN
922 : WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
923 0 : iteration, nmo_converged, max_norm, min_norm, t2 - t1
924 : END IF
925 536 : t1 = m_walltime()
926 :
927 : END DO ! iteration
928 :
929 64 : DEALLOCATE (ritz_coeff)
930 64 : DEALLOCATE (iconv)
931 64 : DEALLOCATE (inotconv)
932 64 : DEALLOCATE (vnorm)
933 :
934 64 : CALL timestop(handle)
935 :
936 192 : END SUBROUTINE generate_extended_space_sparse
937 :
938 : ! **************************************************************************************************
939 :
940 : ! **************************************************************************************************
941 : !> \brief ...
942 : !> \param s_block ...
943 : !> \param h_block ...
944 : !> \param v_block ...
945 : !> \param w_block ...
946 : !> \param evals ...
947 : !> \param ndim ...
948 : ! **************************************************************************************************
949 148 : SUBROUTINE reduce_extended_space(s_block, h_block, v_block, w_block, evals, ndim)
950 :
951 : TYPE(cp_fm_type), INTENT(IN) :: s_block, h_block, v_block, w_block
952 : REAL(dp), DIMENSION(:) :: evals
953 : INTEGER :: ndim
954 :
955 : CHARACTER(len=*), PARAMETER :: routineN = 'reduce_extended_space'
956 :
957 : INTEGER :: handle, info
958 :
959 148 : CALL timeset(routineN, handle)
960 :
961 148 : CALL cp_fm_to_fm(s_block, w_block)
962 148 : CALL cp_fm_cholesky_decompose(s_block, info_out=info)
963 148 : IF (info == 0) THEN
964 148 : CALL cp_fm_triangular_invert(s_block)
965 148 : CALL cp_fm_cholesky_restore(H_block, ndim, S_block, w_block, "MULTIPLY", pos="RIGHT")
966 148 : CALL cp_fm_cholesky_restore(w_block, ndim, S_block, H_block, "MULTIPLY", pos="LEFT", transa="T")
967 148 : CALL choose_eigv_solver(H_block, w_block, evals)
968 148 : CALL cp_fm_cholesky_restore(w_block, ndim, S_block, v_block, "MULTIPLY")
969 : ELSE
970 : ! S^(-1/2)
971 0 : CALL cp_fm_power(w_block, s_block, -0.5_dp, 1.0E-5_dp, info)
972 0 : CALL cp_fm_to_fm(w_block, s_block)
973 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, H_block, s_block, 0.0_dp, w_block)
974 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, H_block)
975 0 : CALL choose_eigv_solver(H_block, w_block, evals)
976 0 : CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, v_block)
977 : END IF
978 :
979 148 : CALL timestop(handle)
980 :
981 148 : END SUBROUTINE reduce_extended_space
982 :
983 : END MODULE qs_scf_block_davidson
|
