Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 itrative
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 30 : 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 30 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
95 30 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
96 : LOGICAL :: converged, do_apply_preconditioner
97 : REAL(dp) :: lambda, max_norm, min_norm, t1, t2
98 30 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm
99 30 : 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 30 : CALL timeset(routineN, handle)
108 :
109 30 : NULLIFY (mo_coeff, mo_coeff_b, eigenvalues)
110 :
111 30 : do_apply_preconditioner = .FALSE.
112 30 : 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 30 : nao=nao, nmo=nmo, homo=homo)
115 30 : IF (do_apply_preconditioner) THEN
116 28 : max_iter = bdav_env%max_iter
117 : ELSE
118 : max_iter = 1
119 : END IF
120 :
121 30 : NULLIFY (c_z, c_pz)
122 30 : NULLIFY (evals, eig_not_conv)
123 30 : t1 = m_walltime()
124 30 : 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 90 : ALLOCATE (iconv(nmo))
130 60 : ALLOCATE (inotconv(nmo))
131 90 : ALLOCATE (ritz_coeff(nmo))
132 60 : ALLOCATE (vnorm(nmo))
133 :
134 30 : converged = .FALSE.
135 82 : DO iter = 1, max_iter
136 :
137 : ! compute Ritz values
138 1998 : ritz_coeff = 0.0_dp
139 54 : CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc")
140 54 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo)
141 54 : CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc")
142 54 : 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 54 : para_env=mo_coeff%matrix_struct%para_env)
147 54 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
148 54 : CALL cp_fm_struct_release(fm_struct_tmp)
149 :
150 54 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp)
151 54 : CALL cp_fm_get_diag(m_tmp, ritz_coeff)
152 54 : CALL cp_fm_release(m_tmp)
153 :
154 : ! Check for converged eigenvectors
155 54 : c_z => bdav_env%matrix_z
156 54 : c_pz => bdav_env%matrix_pz
157 54 : CALL cp_fm_to_fm(m_sc, c_z)
158 54 : CALL cp_fm_column_scale(c_z, ritz_coeff)
159 54 : CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc)
160 54 : CALL cp_fm_vectorsnorm(c_z, vnorm)
161 :
162 54 : nmo_converged = 0
163 54 : nmo_not_converged = 0
164 54 : max_norm = 0.0_dp
165 54 : min_norm = 1.e10_dp
166 1998 : DO imo = 1, nmo
167 1944 : max_norm = MAX(max_norm, vnorm(imo))
168 1998 : min_norm = MIN(min_norm, vnorm(imo))
169 : END DO
170 1998 : iconv = 0
171 1998 : inotconv = 0
172 1998 : DO imo = 1, nmo
173 1998 : IF (vnorm(imo) <= bdav_env%eps_iter) THEN
174 76 : nmo_converged = nmo_converged + 1
175 76 : iconv(nmo_converged) = imo
176 : ELSE
177 1868 : nmo_not_converged = nmo_not_converged + 1
178 1868 : inotconv(nmo_not_converged) = imo
179 : END IF
180 : END DO
181 :
182 54 : IF (nmo_converged > 0) THEN
183 24 : ALLOCATE (iconv_set(nmo_converged, 2))
184 24 : ALLOCATE (inotconv_set(nmo_not_converged, 2))
185 8 : i_last = iconv(1)
186 8 : nset = 0
187 84 : DO j = 1, nmo_converged
188 76 : imo = iconv(j)
189 :
190 84 : IF (imo == i_last + 1) THEN
191 60 : i_last = imo
192 60 : iconv_set(nset, 2) = imo
193 : ELSE
194 16 : i_last = imo
195 16 : nset = nset + 1
196 16 : iconv_set(nset, 1) = imo
197 16 : iconv_set(nset, 2) = imo
198 : END IF
199 : END DO
200 8 : nset_conv = nset
201 :
202 8 : i_last = inotconv(1)
203 8 : nset = 0
204 220 : DO j = 1, nmo_not_converged
205 212 : imo = inotconv(j)
206 :
207 220 : IF (imo == i_last + 1) THEN
208 192 : i_last = imo
209 192 : inotconv_set(nset, 2) = imo
210 : ELSE
211 20 : i_last = imo
212 20 : nset = nset + 1
213 20 : inotconv_set(nset, 1) = imo
214 20 : inotconv_set(nset, 2) = imo
215 : END IF
216 : END DO
217 8 : nset_not_conv = nset
218 8 : CALL cp_fm_release(m_sc)
219 8 : CALL cp_fm_release(m_hc)
220 8 : NULLIFY (c_z, c_pz)
221 : END IF
222 :
223 54 : IF (REAL(nmo_converged, dp)/REAL(nmo, dp) > bdav_env%conv_percent) THEN
224 2 : converged = .TRUE.
225 2 : DEALLOCATE (iconv_set)
226 2 : DEALLOCATE (inotconv_set)
227 2 : t2 = m_walltime()
228 2 : 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 52 : 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 6 : para_env=mo_coeff%matrix_struct%para_env)
243 : !allocate h_fm
244 6 : CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp")
245 : !allocate s_fm
246 6 : CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp")
247 : !copy matrix_h in h_fm
248 6 : CALL copy_dbcsr_to_fm(matrix_h, h_fm)
249 6 : 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 6 : CALL copy_dbcsr_to_fm(matrix_s, s_fm)
254 :
255 : !allocate c_out
256 6 : CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp")
257 : ! set c_out to zero
258 6 : CALL cp_fm_set_all(c_out, 0.0_dp)
259 6 : 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 6 : para_env=mo_coeff%matrix_struct%para_env)
265 6 : CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv")
266 6 : CALL cp_fm_set_all(c_conv, 0.0_dp)
267 : !allocate m_tmp
268 6 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc")
269 6 : 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 52 : para_env=mo_coeff%matrix_struct%para_env)
276 52 : CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv")
277 52 : CALL cp_fm_set_all(c_notconv, 0.0_dp)
278 52 : IF (nmo_converged > 0) THEN
279 6 : CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc")
280 6 : CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc")
281 : !allocate c_z
282 6 : ALLOCATE (c_z, c_pz)
283 6 : CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z")
284 6 : CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz")
285 6 : CALL cp_fm_set_all(c_z, 0.0_dp)
286 :
287 : ! sum contributions to c_out
288 6 : jj = 1
289 16 : DO j = 1, nset_conv
290 10 : i_first = iconv_set(j, 1)
291 10 : i_last = iconv_set(j, 2)
292 10 : n = i_last - i_first + 1
293 10 : CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj)
294 16 : jj = jj + n
295 : END DO
296 6 : CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp)
297 6 : 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 6 : lambda = 100.0_dp*ABS(eigenvalues(homo))
301 6 : CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm)
302 6 : CALL cp_fm_release(m_tmp)
303 6 : CALL cp_fm_release(h_fm)
304 :
305 : END IF
306 :
307 : !allocate m_tmp
308 52 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm")
309 52 : CALL cp_fm_struct_release(fm_struct_tmp)
310 52 : IF (nmo_converged > 0) THEN
311 18 : ALLOCATE (eig_not_conv(nmo_not_converged))
312 6 : jj = 1
313 20 : DO j = 1, nset_not_conv
314 14 : i_first = inotconv_set(j, 1)
315 14 : i_last = inotconv_set(j, 2)
316 14 : n = i_last - i_first + 1
317 14 : CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj)
318 194 : eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
319 20 : jj = jj + n
320 : END DO
321 6 : CALL parallel_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc)
322 6 : 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 6 : CALL cp_fm_to_fm(m_sc, m_tmp)
325 6 : CALL cp_fm_column_scale(m_tmp, eig_not_conv)
326 6 : CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc)
327 6 : DEALLOCATE (eig_not_conv)
328 6 : CALL cp_fm_to_fm(m_tmp, c_z)
329 : ELSE
330 46 : CALL cp_fm_to_fm(mo_coeff, c_notconv)
331 : END IF
332 :
333 : !preconditioner
334 52 : IF (do_apply_preconditioner) THEN
335 50 : IF (preconditioner%in_use /= 0) THEN
336 50 : 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 2 : CALL cp_fm_to_fm(c_z, c_pz)
342 : END IF
343 52 : 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 52 : para_env=mo_coeff%matrix_struct%para_env)
348 :
349 52 : CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm")
350 52 : CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm")
351 52 : CALL cp_fm_struct_release(fm_struct_tmp)
352 :
353 52 : nmat = nmo_not_converged
354 52 : 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 52 : para_env=mo_coeff%matrix_struct%para_env)
358 :
359 52 : CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
360 52 : CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
361 52 : CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
362 52 : CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
363 156 : ALLOCATE (evals(nmat2))
364 :
365 52 : CALL cp_fm_struct_release(fm_struct_tmp)
366 :
367 : ! compute CSC
368 52 : CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
369 :
370 : ! compute CHC
371 52 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp)
372 52 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1)
373 :
374 : ! compute ZSC
375 52 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
376 52 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
377 52 : CALL cp_fm_transpose(m_tmp, mt_tmp)
378 52 : CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
379 : ! compute ZHC
380 52 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
381 52 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
382 52 : CALL cp_fm_transpose(m_tmp, mt_tmp)
383 52 : CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
384 :
385 52 : CALL cp_fm_release(mt_tmp)
386 :
387 : ! reuse m_sc and m_hc to computr HZ and SZ
388 52 : IF (nmo_converged > 0) THEN
389 6 : CALL parallel_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc)
390 6 : CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc)
391 :
392 6 : CALL cp_fm_release(c_out)
393 6 : CALL cp_fm_release(c_conv)
394 6 : CALL cp_fm_release(s_fm)
395 : ELSE
396 46 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo)
397 46 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo)
398 : END IF
399 :
400 : ! compute ZSZ
401 52 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
402 52 : CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
403 : ! compute ZHZ
404 52 : CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
405 52 : CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
406 :
407 52 : CALL cp_fm_release(m_sc)
408 :
409 : ! solution of the reduced eigenvalues problem
410 52 : CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
411 :
412 : ! extract egenvectors
413 52 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1)
414 52 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc)
415 52 : CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1)
416 52 : CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc)
417 :
418 52 : CALL cp_fm_release(m_tmp)
419 :
420 52 : CALL cp_fm_release(c_notconv)
421 52 : CALL cp_fm_release(s_block)
422 52 : CALL cp_fm_release(h_block)
423 52 : CALL cp_fm_release(w_block)
424 52 : CALL cp_fm_release(v_block)
425 :
426 52 : IF (nmo_converged > 0) THEN
427 6 : CALL cp_fm_release(c_z)
428 6 : CALL cp_fm_release(c_pz)
429 6 : DEALLOCATE (c_z, c_pz)
430 6 : jj = 1
431 20 : DO j = 1, nset_not_conv
432 14 : i_first = inotconv_set(j, 1)
433 14 : i_last = inotconv_set(j, 2)
434 14 : n = i_last - i_first + 1
435 14 : CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first)
436 388 : eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
437 20 : jj = jj + n
438 : END DO
439 6 : DEALLOCATE (iconv_set)
440 6 : DEALLOCATE (inotconv_set)
441 : ELSE
442 46 : CALL cp_fm_to_fm(m_hc, mo_coeff)
443 3404 : eigenvalues(1:nmo) = evals(1:nmo)
444 : END IF
445 52 : DEALLOCATE (evals)
446 :
447 52 : CALL cp_fm_release(m_hc)
448 :
449 52 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
450 :
451 52 : t2 = m_walltime()
452 52 : 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 450 : t1 = m_walltime()
457 :
458 : END DO ! iter
459 :
460 30 : DEALLOCATE (iconv)
461 30 : DEALLOCATE (inotconv)
462 30 : DEALLOCATE (ritz_coeff)
463 30 : DEALLOCATE (vnorm)
464 :
465 30 : CALL timestop(handle)
466 90 : 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 3256 : 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 3256 : 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 3256 : iconv = 0
598 3256 : 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 116 : 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 116 : CALL timeset(routineN, handle)
960 :
961 116 : CALL cp_fm_to_fm(s_block, w_block)
962 116 : CALL cp_fm_cholesky_decompose(s_block, info_out=info)
963 116 : IF (info == 0) THEN
964 116 : CALL cp_fm_triangular_invert(s_block)
965 116 : CALL cp_fm_cholesky_restore(H_block, ndim, S_block, w_block, "MULTIPLY", pos="RIGHT")
966 116 : CALL cp_fm_cholesky_restore(w_block, ndim, S_block, H_block, "MULTIPLY", pos="LEFT", transa="T")
967 116 : CALL choose_eigv_solver(H_block, w_block, evals)
968 116 : 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 116 : CALL timestop(handle)
980 :
981 116 : END SUBROUTINE reduce_extended_space
982 :
983 : END MODULE qs_scf_block_davidson
|
