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