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 Apply the direct inversion in the iterative subspace (DIIS) of Pulay
10 : !> in the framework of an SCF iteration for convergence acceleration
11 : !> \par Literature
12 : !> - P. Pulay, Chem. Phys. Lett. 73, 393 (1980)
13 : !> - P. Pulay, J. Comput. Chem. 3, 556 (1982)
14 : !> \par History
15 : !> - Changed to BLACS matrix usage (08.06.2001,MK)
16 : !> - rewritten to include LSD (1st attempt) (01.2003, Joost VandeVondele)
17 : !> - DIIS for ROKS (05.04.06,MK)
18 : !> - DIIS for k-points (04.2023, Augustin Bussy)
19 : !> \author Matthias Krack (28.06.2000)
20 : ! **************************************************************************************************
21 : MODULE qs_diis
22 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm,&
23 : cp_cfm_trace
24 : USE cp_cfm_types, ONLY: cp_cfm_create,&
25 : cp_cfm_get_info,&
26 : cp_cfm_to_cfm,&
27 : cp_cfm_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
30 : dbcsr_set, dbcsr_transposed, dbcsr_type
31 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
32 : dbcsr_maxabs
33 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
35 : cp_fm_scale_and_add,&
36 : cp_fm_symm,&
37 : cp_fm_trace
38 : USE cp_fm_struct, ONLY: cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_get_info,&
41 : cp_fm_maxabsval,&
42 : cp_fm_set_all,&
43 : cp_fm_to_fm,&
44 : cp_fm_type
45 : USE cp_log_handling, ONLY: cp_get_default_logger,&
46 : cp_logger_type,&
47 : cp_to_string
48 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
49 : cp_print_key_unit_nr
50 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
51 : USE input_section_types, ONLY: section_vals_type
52 : USE kinds, ONLY: default_string_length,&
53 : dp
54 : USE local_gemm_api, ONLY: local_gemm_ctxt_type
55 : USE mathlib, ONLY: diag_complex,&
56 : diamat_all
57 : USE message_passing, ONLY: mp_para_env_type
58 : USE parallel_gemm_api, ONLY: parallel_gemm
59 : USE qs_diis_types, ONLY: qs_diis_buffer_type,&
60 : qs_diis_buffer_type_kp,&
61 : qs_diis_buffer_type_sparse
62 : USE qs_environment_types, ONLY: get_qs_env,&
63 : qs_environment_type
64 : USE qs_mo_types, ONLY: get_mo_set,&
65 : mo_set_type
66 : USE string_utilities, ONLY: compress
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_diis'
74 :
75 : ! Public subroutines
76 :
77 : PUBLIC :: qs_diis_b_clear, &
78 : qs_diis_b_create, &
79 : qs_diis_b_step
80 : PUBLIC :: qs_diis_b_clear_sparse, &
81 : qs_diis_b_create_sparse, &
82 : qs_diis_b_step_4lscf
83 : PUBLIC :: qs_diis_b_clear_kp, &
84 : qs_diis_b_check_i_alloc_kp, &
85 : qs_diis_b_create_kp, &
86 : qs_diis_b_step_kp, &
87 : qs_diis_b_calc_err_kp, &
88 : qs_diis_b_info_kp
89 :
90 : CONTAINS
91 :
92 : ! **************************************************************************************************
93 : !> \brief Allocates an SCF DIIS buffer
94 : !> \param diis_buffer the buffer to create
95 : !> \param nbuffer ...
96 : !> \par History
97 : !> 02.2003 created [fawzi]
98 : !> \author fawzi
99 : ! **************************************************************************************************
100 4608 : SUBROUTINE qs_diis_b_create(diis_buffer, nbuffer)
101 :
102 : TYPE(qs_diis_buffer_type), INTENT(OUT) :: diis_buffer
103 : INTEGER, INTENT(in) :: nbuffer
104 :
105 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_diis_b_create'
106 :
107 : INTEGER :: handle
108 :
109 : ! -------------------------------------------------------------------------
110 :
111 4608 : CALL timeset(routineN, handle)
112 :
113 4608 : NULLIFY (diis_buffer%b_matrix)
114 4608 : NULLIFY (diis_buffer%error)
115 4608 : NULLIFY (diis_buffer%param)
116 4608 : diis_buffer%nbuffer = nbuffer
117 4608 : diis_buffer%ncall = 0
118 :
119 4608 : CALL timestop(handle)
120 :
121 4608 : END SUBROUTINE qs_diis_b_create
122 :
123 : ! **************************************************************************************************
124 : !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
125 : !> variables and with a buffer size of nbuffer.
126 : !> \param diis_buffer the buffer to initialize
127 : !> \param matrix_struct the structure for the matrix of the buffer
128 : !> \param nspin ...
129 : !> \param scf_section ...
130 : !> \par History
131 : !> - Creation (07.05.2001, Matthias Krack)
132 : !> - Changed to BLACS matrix usage (08.06.2001,MK)
133 : !> - DIIS for ROKS (05.04.06,MK)
134 : !> \author Matthias Krack
135 : !> \note
136 : !> check to allocate matrixes only when needed, using a linked list?
137 : ! **************************************************************************************************
138 99409 : SUBROUTINE qs_diis_b_check_i_alloc(diis_buffer, matrix_struct, nspin, &
139 : scf_section)
140 :
141 : TYPE(qs_diis_buffer_type), INTENT(INOUT) :: diis_buffer
142 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
143 : INTEGER, INTENT(IN) :: nspin
144 : TYPE(section_vals_type), POINTER :: scf_section
145 :
146 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc'
147 :
148 : INTEGER :: handle, ibuffer, ispin, nbuffer, &
149 : output_unit
150 : TYPE(cp_logger_type), POINTER :: logger
151 :
152 : ! -------------------------------------------------------------------------
153 :
154 99409 : CALL timeset(routineN, handle)
155 :
156 99409 : logger => cp_get_default_logger()
157 :
158 99409 : nbuffer = diis_buffer%nbuffer
159 :
160 99409 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
161 35430 : ALLOCATE (diis_buffer%error(nbuffer, nspin))
162 :
163 7792 : DO ispin = 1, nspin
164 24840 : DO ibuffer = 1, nbuffer
165 : CALL cp_fm_create(diis_buffer%error(ibuffer, ispin), &
166 : name="qs_diis_b%error("// &
167 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
168 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
169 21310 : matrix_struct=matrix_struct)
170 : END DO
171 : END DO
172 : END IF
173 :
174 99409 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
175 35430 : ALLOCATE (diis_buffer%param(nbuffer, nspin))
176 :
177 7792 : DO ispin = 1, nspin
178 24840 : DO ibuffer = 1, nbuffer
179 : CALL cp_fm_create(diis_buffer%param(ibuffer, ispin), &
180 : name="qs_diis_b%param("// &
181 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
182 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
183 21310 : matrix_struct=matrix_struct)
184 : END DO
185 : END DO
186 : END IF
187 :
188 99409 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
189 14120 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
190 109430 : diis_buffer%b_matrix = 0.0_dp
191 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
192 3530 : extension=".scfLog")
193 3530 : IF (output_unit > 0) THEN
194 : WRITE (UNIT=output_unit, FMT="(/,T9,A)") &
195 20 : "DIIS | The SCF DIIS buffer was allocated and initialized"
196 : END IF
197 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
198 3530 : "PRINT%DIIS_INFO")
199 : END IF
200 :
201 99409 : CALL timestop(handle)
202 :
203 99409 : END SUBROUTINE qs_diis_b_check_i_alloc
204 :
205 : ! **************************************************************************************************
206 : !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step.
207 : !> \param diis_buffer ...
208 : !> \param mo_array ...
209 : !> \param kc ...
210 : !> \param sc ...
211 : !> \param delta ...
212 : !> \param error_max ...
213 : !> \param diis_step ...
214 : !> \param eps_diis ...
215 : !> \param nmixing ...
216 : !> \param s_matrix ...
217 : !> \param scf_section ...
218 : !> \param roks ...
219 : !> \par History
220 : !> - Creation (07.05.2001, Matthias Krack)
221 : !> - Changed to BLACS matrix usage (08.06.2001, MK)
222 : !> - 03.2003 rewamped [fawzi]
223 : !> - Adapted for high-spin ROKS (08.04.06,MK)
224 : !> \author Matthias Krack
225 : ! **************************************************************************************************
226 99409 : SUBROUTINE qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, &
227 : diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
228 :
229 : TYPE(qs_diis_buffer_type), POINTER :: diis_buffer
230 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
231 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: kc
232 : TYPE(cp_fm_type), INTENT(IN) :: sc
233 : REAL(KIND=dp), INTENT(IN) :: delta
234 : REAL(KIND=dp), INTENT(OUT) :: error_max
235 : LOGICAL, INTENT(OUT) :: diis_step
236 : REAL(KIND=dp), INTENT(IN) :: eps_diis
237 : INTEGER, INTENT(IN), OPTIONAL :: nmixing
238 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
239 : POINTER :: s_matrix
240 : TYPE(section_vals_type), POINTER :: scf_section
241 : LOGICAL, INTENT(IN), OPTIONAL :: roks
242 :
243 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step'
244 : REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
245 :
246 : CHARACTER(LEN=2*default_string_length) :: message
247 : INTEGER :: handle, homo, ib, imo, ispin, jb, &
248 : my_nmixing, nao, nb, nb1, nmo, nspin, &
249 : output_unit
250 : LOGICAL :: eigenvectors_discarded, my_roks
251 : REAL(KIND=dp) :: maxocc, tmp
252 99409 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev, occ, traces
253 99409 : REAL(KIND=dp), DIMENSION(:), POINTER :: occa, occb
254 99409 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
255 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
256 : TYPE(cp_fm_type), POINTER :: c, new_errors, parameters
257 : TYPE(cp_logger_type), POINTER :: logger
258 :
259 : ! -------------------------------------------------------------------------
260 :
261 99409 : CALL timeset(routineN, handle)
262 :
263 99409 : nspin = SIZE(mo_array)
264 99409 : diis_step = .FALSE.
265 :
266 99409 : IF (PRESENT(roks)) THEN
267 1046 : my_roks = .TRUE.
268 1046 : nspin = 1
269 : ELSE
270 : my_roks = .FALSE.
271 : END IF
272 :
273 99409 : my_nmixing = 2
274 99409 : IF (PRESENT(nmixing)) my_nmixing = nmixing
275 :
276 99409 : NULLIFY (c, new_errors, parameters, matrix_struct, a, b, occa, occb)
277 99409 : logger => cp_get_default_logger()
278 :
279 : ! Quick return, if no DIIS is requested
280 :
281 99409 : IF (diis_buffer%nbuffer < 1) THEN
282 0 : CALL timestop(handle)
283 0 : RETURN
284 : END IF
285 :
286 : CALL cp_fm_get_info(kc(1), &
287 99409 : matrix_struct=matrix_struct)
288 : CALL qs_diis_b_check_i_alloc(diis_buffer, &
289 : matrix_struct=matrix_struct, &
290 : nspin=nspin, &
291 99409 : scf_section=scf_section)
292 :
293 99409 : error_max = 0.0_dp
294 :
295 99409 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
296 99409 : diis_buffer%ncall = diis_buffer%ncall + 1
297 99409 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
298 :
299 99409 : IF (PRESENT(s_matrix)) THEN
300 83667 : IF (.NOT. ALLOCATED(diis_buffer%overlap_copy_plan)) ALLOCATE (diis_buffer%overlap_copy_plan)
301 : END IF
302 :
303 214976 : DO ispin = 1, nspin
304 :
305 : CALL get_mo_set(mo_set=mo_array(ispin), &
306 : nao=nao, &
307 : nmo=nmo, &
308 : homo=homo, &
309 : mo_coeff=c, &
310 : occupation_numbers=occa, &
311 115567 : maxocc=maxocc)
312 :
313 115567 : new_errors => diis_buffer%error(ib, ispin)
314 115567 : parameters => diis_buffer%param(ib, ispin)
315 :
316 : ! Copy the Kohn-Sham matrix K to the DIIS buffer
317 :
318 115567 : CALL cp_fm_to_fm(kc(ispin), parameters)
319 :
320 115567 : IF (my_roks) THEN
321 :
322 3138 : ALLOCATE (occ(nmo))
323 :
324 : CALL get_mo_set(mo_set=mo_array(2), &
325 1046 : occupation_numbers=occb)
326 :
327 16114 : DO imo = 1, nmo
328 16114 : occ(imo) = SQRT(occa(imo) + occb(imo))
329 : END DO
330 :
331 1046 : CALL cp_fm_to_fm(c, sc)
332 1046 : CALL cp_fm_column_scale(sc, occ(1:homo))
333 :
334 : ! KC <- K*C
335 1046 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, parameters, sc, 0.0_dp, kc(ispin))
336 :
337 1046 : IF (PRESENT(s_matrix)) THEN
338 586 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors, diis_buffer%overlap_copy_plan, check_plan=.TRUE.)
339 : ! SC <- S*C
340 586 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, new_errors, c, 0.0_dp, sc)
341 586 : CALL cp_fm_column_scale(sc, occ(1:homo))
342 : END IF
343 :
344 : ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
345 : ! or for an orthogonal basis
346 : ! new_errors <- KC*C^T - C*(KC)^T = K*P - P*K with S = I
347 1046 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
348 1046 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
349 :
350 1046 : DEALLOCATE (occ)
351 :
352 : ELSE
353 :
354 : ! KC <- K*C
355 114521 : CALL cp_fm_symm("L", "U", nao, homo, maxocc, parameters, c, 0.0_dp, kc(ispin))
356 :
357 114521 : IF (PRESENT(s_matrix)) THEN
358 : ! I guess that this copy can be avoided for LSD
359 98545 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors, diis_buffer%overlap_copy_plan, check_plan=.TRUE.)
360 : ! sc <- S*C
361 98545 : CALL cp_fm_symm("L", "U", nao, homo, 2.0_dp, new_errors, c, 0.0_dp, sc)
362 : ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
363 98545 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
364 98545 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
365 : ELSE
366 : ! new_errors <- KC*(C)^T - C*(KC)^T = K*P - P*K
367 15976 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, c, kc(ispin), 0.0_dp, new_errors)
368 15976 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), c, -1.0_dp, new_errors)
369 : END IF
370 :
371 : END IF
372 :
373 115567 : CALL cp_fm_maxabsval(new_errors, tmp)
374 330543 : error_max = MAX(error_max, tmp)
375 :
376 : END DO
377 :
378 : ! Check, if a DIIS step is appropriate
379 :
380 99409 : diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
381 :
382 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
383 99409 : extension=".scfLog")
384 99409 : IF (output_unit > 0) THEN
385 : WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
386 254 : "DIIS | Current SCF DIIS buffer size: ", nb, &
387 254 : "DIIS | Maximum SCF DIIS error vector element:", error_max, &
388 254 : "DIIS | Current SCF convergence: ", delta, &
389 508 : "DIIS | Threshold value for a DIIS step: ", eps_diis
390 254 : IF (error_max < eps_diis) THEN
391 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
392 102 : "DIIS | => The SCF DIIS buffer will be updated"
393 : ELSE
394 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
395 152 : "DIIS | => No update of the SCF DIIS buffer"
396 : END IF
397 254 : IF (diis_step .AND. (error_max < eps_diis)) THEN
398 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
399 62 : "DIIS | => A SCF DIIS step will be performed"
400 : ELSE
401 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
402 192 : "DIIS | => No SCF DIIS step will be performed"
403 : END IF
404 : END IF
405 :
406 : ! Update the SCF DIIS buffer
407 :
408 99409 : IF (error_max < eps_diis) THEN
409 :
410 86139 : b => diis_buffer%b_matrix
411 :
412 258417 : ALLOCATE (traces(nb))
413 365867 : b(1:nb, ib) = 0.0_dp
414 : ! Batch history traces, but avoid starting OpenMP for a single entry.
415 186584 : DO ispin = 1, nspin
416 100445 : IF (nb == 1) THEN
417 11552 : CALL cp_fm_trace(diis_buffer%error(1, ispin), diis_buffer%error(ib, ispin), traces(1))
418 : ELSE
419 88893 : CALL cp_fm_trace(diis_buffer%error(1:nb, ispin), diis_buffer%error(ib, ispin), traces)
420 : END IF
421 515112 : b(1:nb, ib) = b(1:nb, ib) + traces(:)
422 : END DO
423 86139 : DEALLOCATE (traces)
424 :
425 365867 : DO jb = 1, nb
426 365867 : b(ib, jb) = b(jb, ib)
427 : END DO
428 :
429 : ELSE
430 :
431 13270 : diis_step = .FALSE.
432 :
433 : END IF
434 :
435 : ! Perform DIIS step
436 :
437 99409 : IF (diis_step) THEN
438 :
439 64257 : nb1 = nb + 1
440 :
441 257028 : ALLOCATE (a(nb1, nb1))
442 192771 : ALLOCATE (b(nb1, nb1))
443 192771 : ALLOCATE (ev(nb1))
444 :
445 : ! Set up the linear DIIS equation system
446 :
447 2267217 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
448 :
449 296255 : b(1:nb, nb1) = -1.0_dp
450 296255 : b(nb1, 1:nb) = -1.0_dp
451 64257 : b(nb1, nb1) = 0.0_dp
452 :
453 : ! Solve the linear DIIS equation system
454 :
455 360512 : ev(1:nb1) = 0.0_dp
456 64257 : CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
457 :
458 3452237 : a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
459 :
460 64257 : eigenvectors_discarded = .FALSE.
461 :
462 360512 : DO jb = 1, nb1
463 360512 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
464 45076 : IF (output_unit > 0) THEN
465 5 : IF (.NOT. eigenvectors_discarded) THEN
466 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
467 5 : "DIIS | Checking eigenvalues of the DIIS error matrix"
468 : END IF
469 : WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
470 5 : "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
471 10 : "threshold ", eigenvalue_threshold
472 5 : CALL compress(message)
473 5 : WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
474 5 : eigenvectors_discarded = .TRUE.
475 : END IF
476 263932 : a(1:nb1, jb) = 0.0_dp
477 : ELSE
478 1430058 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
479 : END IF
480 : END DO
481 :
482 64257 : IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
483 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
484 5 : "DIIS | The corresponding eigenvectors were discarded"
485 : END IF
486 :
487 1693990 : ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
488 :
489 : ! Update Kohn-Sham matrix
490 :
491 139334 : DO ispin = 1, nspin
492 75077 : CALL cp_fm_set_all(kc(ispin), 0.0_dp)
493 411890 : DO jb = 1, nb
494 272556 : parameters => diis_buffer%param(jb, ispin)
495 347633 : CALL cp_fm_scale_and_add(1.0_dp, kc(ispin), -ev(jb), parameters)
496 : END DO
497 : END DO
498 :
499 64257 : DEALLOCATE (a)
500 64257 : DEALLOCATE (b)
501 64257 : DEALLOCATE (ev)
502 :
503 : ELSE
504 :
505 75642 : DO ispin = 1, nspin
506 40490 : parameters => diis_buffer%param(ib, ispin)
507 75642 : CALL cp_fm_to_fm(parameters, kc(ispin))
508 : END DO
509 :
510 : END IF
511 :
512 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
513 99409 : "PRINT%DIIS_INFO")
514 :
515 99409 : CALL timestop(handle)
516 :
517 198818 : END SUBROUTINE qs_diis_b_step
518 :
519 : ! **************************************************************************************************
520 : !> \brief clears the buffer
521 : !> \param diis_buffer the buffer to clear
522 : !> \par History
523 : !> 02.2003 created [fawzi]
524 : !> \author fawzi
525 : ! **************************************************************************************************
526 16586 : PURE SUBROUTINE qs_diis_b_clear(diis_buffer)
527 :
528 : TYPE(qs_diis_buffer_type), INTENT(INOUT) :: diis_buffer
529 :
530 16586 : diis_buffer%ncall = 0
531 :
532 16586 : END SUBROUTINE qs_diis_b_clear
533 :
534 : ! **************************************************************************************************
535 : !> \brief Update the SCF DIIS buffer in linear scaling SCF (LS-SCF),
536 : !> and if appropriate does a diis step.
537 : !> \param diis_buffer ...
538 : !> \param qs_env ...
539 : !> \param ls_scf_env ...
540 : !> \param unit_nr ...
541 : !> \param iscf ...
542 : !> \param diis_step ...
543 : !> \param eps_diis ...
544 : !> \param nmixing ...
545 : !> \param s_matrix ...
546 : !> \param threshold ...
547 : !> \par History
548 : !> - Adapted for LS-SCF (10-11-14) from qs_diis_b_step
549 : !> \author Fredy W. Aquino
550 : ! **************************************************************************************************
551 :
552 18 : SUBROUTINE qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, &
553 : diis_step, eps_diis, nmixing, s_matrix, threshold)
554 : ! Note.- Input: ls_scf_env%matrix_p(ispin) , Density Matrix
555 : ! matrix_ks (from qs_env) , Kohn-Sham Matrix (IN/OUT)
556 :
557 : TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer
558 : TYPE(qs_environment_type), POINTER :: qs_env
559 : TYPE(ls_scf_env_type) :: ls_scf_env
560 : INTEGER, INTENT(IN) :: unit_nr, iscf
561 : LOGICAL, INTENT(OUT) :: diis_step
562 : REAL(KIND=dp), INTENT(IN) :: eps_diis
563 : INTEGER, INTENT(IN), OPTIONAL :: nmixing
564 : TYPE(dbcsr_type), OPTIONAL :: s_matrix
565 : REAL(KIND=dp), INTENT(IN) :: threshold
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_4lscf'
568 : REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
569 :
570 : INTEGER :: handle, ib, ispin, jb, my_nmixing, nb, &
571 : nb1, nspin
572 : REAL(KIND=dp) :: error_max, tmp
573 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
574 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
575 : TYPE(cp_logger_type), POINTER :: logger
576 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
577 : TYPE(dbcsr_type) :: matrix_KSerr_t, matrix_tmp
578 : TYPE(dbcsr_type), POINTER :: new_errors, old_errors, parameters
579 : TYPE(mp_para_env_type), POINTER :: para_env
580 :
581 18 : CALL timeset(routineN, handle)
582 18 : IF (ls_scf_env%do_pao) THEN
583 0 : CPABORT("LS_SCF%LS_DIIS not compatible with PAO")
584 : END IF
585 18 : nspin = ls_scf_env%nspins
586 18 : diis_step = .FALSE.
587 18 : my_nmixing = 2
588 18 : IF (PRESENT(nmixing)) my_nmixing = nmixing
589 18 : NULLIFY (new_errors, old_errors, parameters, a, b)
590 18 : logger => cp_get_default_logger()
591 : ! Quick return, if no DIIS is requested
592 18 : IF (diis_buffer%nbuffer < 1) THEN
593 0 : CALL timestop(handle)
594 : RETURN
595 : END IF
596 :
597 : ! Getting current Kohn-Sham matrix from qs_env
598 : CALL get_qs_env(qs_env, &
599 : para_env=para_env, &
600 18 : matrix_ks=matrix_ks)
601 : CALL qs_diis_b_check_i_alloc_sparse( &
602 : diis_buffer, &
603 : ls_scf_env, &
604 18 : nspin)
605 18 : error_max = 0.0_dp
606 :
607 18 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
608 18 : diis_buffer%ncall = diis_buffer%ncall + 1
609 18 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
610 : ! Create scratch arrays
611 : CALL dbcsr_create(matrix_tmp, &
612 : template=ls_scf_env%matrix_ks(1), &
613 18 : matrix_type='N')
614 18 : CALL dbcsr_set(matrix_tmp, 0.0_dp) ! reset matrix
615 : CALL dbcsr_create(matrix_KSerr_t, &
616 : template=ls_scf_env%matrix_ks(1), &
617 18 : matrix_type='N')
618 18 : CALL dbcsr_set(matrix_KSerr_t, 0.0_dp) ! reset matrix
619 :
620 46 : DO ispin = 1, nspin ! ------ Loop-ispin----START
621 :
622 28 : new_errors => diis_buffer%error(ib, ispin)%matrix
623 28 : parameters => diis_buffer%param(ib, ispin)%matrix
624 : ! Copy the Kohn-Sham matrix K to the DIIS buffer
625 : CALL dbcsr_copy(parameters, & ! out
626 28 : matrix_ks(ispin)%matrix) ! in
627 :
628 28 : IF (PRESENT(s_matrix)) THEN ! if-s_matrix ---------- START
629 : ! Calculate Kohn-Sham error (non-orthogonal)= K*P*S-(K*P*S)^T
630 : ! matrix_tmp = P*S
631 : CALL dbcsr_multiply("N", "N", &
632 : 1.0_dp, ls_scf_env%matrix_p(ispin), &
633 : s_matrix, &
634 : 0.0_dp, matrix_tmp, &
635 28 : filter_eps=threshold)
636 : ! new_errors= K*P*S
637 : CALL dbcsr_multiply("N", "N", &
638 : 1.0_dp, matrix_ks(ispin)%matrix, &
639 : matrix_tmp, &
640 : 0.0_dp, new_errors, &
641 28 : filter_eps=threshold)
642 : ! matrix_KSerr_t= transpose(K*P*S)
643 : CALL dbcsr_transposed(matrix_KSerr_t, &
644 28 : new_errors)
645 : ! new_errors=K*P*S-transpose(K*P*S)
646 : CALL dbcsr_add(new_errors, &
647 : matrix_KSerr_t, &
648 28 : 1.0_dp, -1.0_dp)
649 : ELSE ! if-s_matrix ---------- MID
650 : ! Calculate Kohn-Sham error (orthogonal)= K*P - P*K
651 : ! new_errors=K*P
652 : CALL dbcsr_multiply("N", "N", &
653 : 1.0_dp, matrix_ks(ispin)%matrix, &
654 : ls_scf_env%matrix_p(ispin), &
655 : 0.0_dp, new_errors, &
656 0 : filter_eps=threshold)
657 : ! matrix_KSerr_t= transpose(K*P)
658 : CALL dbcsr_transposed(matrix_KSerr_t, &
659 0 : new_errors)
660 : ! new_errors=K*P-transpose(K*P)
661 : CALL dbcsr_add(new_errors, &
662 : matrix_KSerr_t, &
663 0 : 1.0_dp, -1.0_dp)
664 : END IF ! if-s_matrix ---------- END
665 :
666 28 : tmp = dbcsr_maxabs(new_errors)
667 46 : error_max = MAX(error_max, tmp)
668 :
669 : END DO ! ------ Loop-ispin----END
670 :
671 : ! Check, if a DIIS step is appropriate
672 :
673 18 : diis_step = (diis_buffer%ncall >= my_nmixing)
674 :
675 18 : IF (unit_nr > 0) THEN
676 : WRITE (unit_nr, '(A29,I3,A3,4(I3,A1))') &
677 9 : "DIIS: (ncall,nbuffer,ib,nb)=(", iscf, ")=(", &
678 18 : diis_buffer%ncall, ",", diis_buffer%nbuffer, ",", ib, ",", nb, ")"
679 : WRITE (unit_nr, '(A57,I3,A3,L1,A1,F10.8,A1,F4.2,A1,L1,A1)') &
680 9 : "DIIS: (diis_step,error_max,eps_diis,error_max<eps_diis)=(", &
681 9 : iscf, ")=(", diis_step, ",", error_max, ",", eps_diis, ",", &
682 18 : (error_max < eps_diis), ")"
683 : WRITE (unit_nr, '(A75)') &
684 9 : "DIIS: diis_step=T : Perform DIIS error_max<eps_diis=T : Update DIIS buffer"
685 : END IF
686 :
687 : ! Update the SCF DIIS buffer
688 18 : IF (error_max < eps_diis) THEN
689 18 : b => diis_buffer%b_matrix
690 66 : DO jb = 1, nb
691 48 : b(jb, ib) = 0.0_dp
692 124 : DO ispin = 1, nspin
693 76 : old_errors => diis_buffer%error(jb, ispin)%matrix
694 76 : new_errors => diis_buffer%error(ib, ispin)%matrix
695 : CALL dbcsr_dot(old_errors, &
696 : new_errors, &
697 76 : tmp) ! out : < f_i | f_j >
698 124 : b(jb, ib) = b(jb, ib) + tmp
699 : END DO ! end-loop-ispin
700 66 : b(ib, jb) = b(jb, ib)
701 : END DO ! end-loop-jb
702 : ELSE
703 0 : diis_step = .FALSE.
704 : END IF
705 :
706 : ! Perform DIIS step
707 18 : IF (diis_step) THEN
708 14 : nb1 = nb + 1
709 56 : ALLOCATE (a(nb1, nb1))
710 42 : ALLOCATE (b(nb1, nb1))
711 42 : ALLOCATE (ev(nb1))
712 : ! Set up the linear DIIS equation system
713 398 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
714 58 : b(1:nb, nb1) = -1.0_dp
715 58 : b(nb1, 1:nb) = -1.0_dp
716 14 : b(nb1, nb1) = 0.0_dp
717 : ! Solve the linear DIIS equation system
718 14 : CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
719 630 : a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
720 72 : DO jb = 1, nb1
721 72 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
722 0 : a(1:nb1, jb) = 0.0_dp
723 : ELSE
724 308 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
725 : END IF
726 : END DO ! end-loop-jb
727 :
728 308 : ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
729 :
730 : ! Update Kohn-Sham matrix
731 14 : IF (iscf >= ls_scf_env%iter_ini_diis) THEN ! if-iscf-to-updateKS------ START
732 :
733 14 : IF (unit_nr > 0) THEN
734 7 : WRITE (unit_nr, '(A40,I3)') 'DIIS: Updating Kohn-Sham matrix at iscf=', iscf
735 : END IF
736 :
737 36 : DO ispin = 1, nspin
738 : CALL dbcsr_set(matrix_ks(ispin)%matrix, & ! reset matrix
739 22 : 0.0_dp)
740 106 : DO jb = 1, nb
741 70 : parameters => diis_buffer%param(jb, ispin)%matrix
742 : CALL dbcsr_add(matrix_ks(ispin)%matrix, parameters, &
743 92 : 1.0_dp, -ev(jb))
744 : END DO ! end-loop-jb
745 : END DO ! end-loop-ispin
746 : END IF ! if-iscf-to-updateKS------ END
747 :
748 14 : DEALLOCATE (a)
749 14 : DEALLOCATE (b)
750 14 : DEALLOCATE (ev)
751 :
752 : ELSE
753 10 : DO ispin = 1, nspin
754 6 : parameters => diis_buffer%param(ib, ispin)%matrix
755 : CALL dbcsr_copy(parameters, & ! out
756 10 : matrix_ks(ispin)%matrix) ! in
757 : END DO ! end-loop-ispin
758 : END IF
759 18 : CALL dbcsr_release(matrix_tmp)
760 18 : CALL dbcsr_release(matrix_KSerr_t)
761 18 : CALL timestop(handle)
762 :
763 18 : END SUBROUTINE qs_diis_b_step_4lscf
764 :
765 : ! **************************************************************************************************
766 : !> \brief Allocate and initialize a DIIS buffer with a buffer size of nbuffer.
767 : !> \param diis_buffer the buffer to initialize
768 : !> \param ls_scf_env ...
769 : !> \param nspin ...
770 : !> \par History
771 : !> - Adapted from qs_diis_b_check_i_alloc for sparse matrices and
772 : !> used in LS-SCF module (ls_scf_main) (10-11-14)
773 : !> \author Fredy W. Aquino
774 : !> \note
775 : !> check to allocate matrices only when needed
776 : ! **************************************************************************************************
777 :
778 18 : SUBROUTINE qs_diis_b_check_i_alloc_sparse(diis_buffer, ls_scf_env, &
779 : nspin)
780 :
781 : TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
782 : TYPE(ls_scf_env_type) :: ls_scf_env
783 : INTEGER, INTENT(IN) :: nspin
784 :
785 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_sparse'
786 :
787 : INTEGER :: handle, ibuffer, ispin, nbuffer
788 : TYPE(cp_logger_type), POINTER :: logger
789 :
790 : ! -------------------------------------------------------------------------
791 :
792 18 : CALL timeset(routineN, handle)
793 :
794 18 : logger => cp_get_default_logger()
795 :
796 18 : nbuffer = diis_buffer%nbuffer
797 :
798 18 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
799 46 : ALLOCATE (diis_buffer%error(nbuffer, nspin))
800 :
801 10 : DO ispin = 1, nspin
802 34 : DO ibuffer = 1, nbuffer
803 24 : ALLOCATE (diis_buffer%error(ibuffer, ispin)%matrix)
804 :
805 : CALL dbcsr_create(diis_buffer%error(ibuffer, ispin)%matrix, &
806 : template=ls_scf_env%matrix_ks(1), &
807 30 : matrix_type='N')
808 : END DO
809 : END DO
810 : END IF
811 :
812 18 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
813 46 : ALLOCATE (diis_buffer%param(nbuffer, nspin))
814 :
815 10 : DO ispin = 1, nspin
816 34 : DO ibuffer = 1, nbuffer
817 24 : ALLOCATE (diis_buffer%param(ibuffer, ispin)%matrix)
818 : CALL dbcsr_create(diis_buffer%param(ibuffer, ispin)%matrix, &
819 : template=ls_scf_env%matrix_ks(1), &
820 30 : matrix_type='N')
821 : END DO
822 : END DO
823 : END IF
824 :
825 18 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
826 16 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
827 :
828 124 : diis_buffer%b_matrix = 0.0_dp
829 : END IF
830 :
831 18 : CALL timestop(handle)
832 :
833 18 : END SUBROUTINE qs_diis_b_check_i_alloc_sparse
834 :
835 : ! **************************************************************************************************
836 : !> \brief clears the DIIS buffer in LS-SCF calculation
837 : !> \param diis_buffer the buffer to clear
838 : !> \par History
839 : !> 10-11-14 created [FA] modified from qs_diis_b_clear
840 : !> \author Fredy W. Aquino
841 : ! **************************************************************************************************
842 :
843 4 : PURE SUBROUTINE qs_diis_b_clear_sparse(diis_buffer)
844 :
845 : TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
846 :
847 4 : diis_buffer%ncall = 0
848 :
849 4 : END SUBROUTINE qs_diis_b_clear_sparse
850 :
851 : ! **************************************************************************************************
852 : !> \brief Allocates an SCF DIIS buffer for LS-SCF calculation
853 : !> \param diis_buffer the buffer to create
854 : !> \param nbuffer ...
855 : !> \par History
856 : !> 10-11-14 created [FA] modified from qs_diis_b_create
857 : !> \author Fredy W. Aquino
858 : ! **************************************************************************************************
859 4 : PURE SUBROUTINE qs_diis_b_create_sparse(diis_buffer, nbuffer)
860 :
861 : TYPE(qs_diis_buffer_type_sparse), INTENT(OUT) :: diis_buffer
862 : INTEGER, INTENT(in) :: nbuffer
863 :
864 : NULLIFY (diis_buffer%b_matrix)
865 : NULLIFY (diis_buffer%error)
866 : NULLIFY (diis_buffer%param)
867 4 : diis_buffer%nbuffer = nbuffer
868 4 : diis_buffer%ncall = 0
869 :
870 4 : END SUBROUTINE qs_diis_b_create_sparse
871 :
872 : ! **************************************************************************************************
873 : !> \brief Allocates an SCF DIIS buffer for k-points
874 : !> \param diis_buffer the buffer to create
875 : !> \param nbuffer ...
876 : ! **************************************************************************************************
877 2560 : SUBROUTINE qs_diis_b_create_kp(diis_buffer, nbuffer)
878 :
879 : TYPE(qs_diis_buffer_type_kp), INTENT(OUT) :: diis_buffer
880 : INTEGER, INTENT(in) :: nbuffer
881 :
882 : NULLIFY (diis_buffer%b_matrix)
883 : NULLIFY (diis_buffer%error)
884 : NULLIFY (diis_buffer%param)
885 : NULLIFY (diis_buffer%smat)
886 2560 : diis_buffer%nbuffer = nbuffer
887 2560 : diis_buffer%ncall = 0
888 :
889 2560 : END SUBROUTINE qs_diis_b_create_kp
890 :
891 : ! **************************************************************************************************
892 : !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
893 : !> variables and with a buffer size of nbuffer, in the k-point case
894 : !> \param diis_buffer the buffer to initialize
895 : !> \param matrix_struct the structure for the matrix of the buffer note: this is in the kp subgroup
896 : !> \param nspin ...
897 : !> \param nkp ...
898 : !> \param scf_section ...
899 : ! **************************************************************************************************
900 27450 : SUBROUTINE qs_diis_b_check_i_alloc_kp(diis_buffer, matrix_struct, nspin, nkp, scf_section)
901 :
902 : TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
903 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
904 : INTEGER, INTENT(IN) :: nspin, nkp
905 : TYPE(section_vals_type), POINTER :: scf_section
906 :
907 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_kp'
908 :
909 : INTEGER :: handle, ibuffer, ikp, ispin, nbuffer, &
910 : output_unit
911 : TYPE(cp_logger_type), POINTER :: logger
912 :
913 : ! -------------------------------------------------------------------------
914 :
915 27450 : CALL timeset(routineN, handle)
916 :
917 27450 : logger => cp_get_default_logger()
918 :
919 27450 : nbuffer = diis_buffer%nbuffer
920 :
921 27450 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
922 35735 : ALLOCATE (diis_buffer%error(nbuffer, nspin, nkp))
923 :
924 5927 : DO ikp = 1, nkp
925 10231 : DO ispin = 1, nspin
926 25375 : DO ibuffer = 1, nbuffer
927 : CALL cp_cfm_create(diis_buffer%error(ibuffer, ispin, ikp), &
928 : name="qs_diis_b%error("// &
929 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
930 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
931 21520 : matrix_struct=matrix_struct)
932 : END DO
933 : END DO
934 : END DO
935 : END IF
936 :
937 27450 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
938 35735 : ALLOCATE (diis_buffer%param(nbuffer, nspin, nkp))
939 :
940 5927 : DO ikp = 1, nkp
941 10231 : DO ispin = 1, nspin
942 25375 : DO ibuffer = 1, nbuffer
943 : CALL cp_cfm_create(diis_buffer%param(ibuffer, ispin, ikp), &
944 : name="qs_diis_b%param("// &
945 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
946 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
947 21520 : matrix_struct=matrix_struct)
948 : END DO
949 : END DO
950 : END DO
951 : END IF
952 :
953 27450 : IF (.NOT. ASSOCIATED(diis_buffer%smat)) THEN
954 10071 : ALLOCATE (diis_buffer%smat(nkp))
955 5927 : DO ikp = 1, nkp
956 : CALL cp_cfm_create(diis_buffer%smat(ikp), &
957 : name="kp_cfm_smat("//TRIM(cp_to_string(ikp))//")", &
958 5927 : matrix_struct=matrix_struct)
959 : END DO
960 : END IF
961 :
962 27450 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
963 8288 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
964 64232 : diis_buffer%b_matrix = 0.0_dp
965 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
966 2072 : extension=".scfLog")
967 2072 : IF (output_unit > 0) THEN
968 : WRITE (UNIT=output_unit, FMT="(/,T9,A)") &
969 0 : "DIIS | The SCF DIIS buffer was allocated and initialized"
970 : END IF
971 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
972 2072 : "PRINT%DIIS_INFO")
973 : END IF
974 :
975 27450 : CALL timestop(handle)
976 :
977 27450 : END SUBROUTINE qs_diis_b_check_i_alloc_kp
978 :
979 : ! **************************************************************************************************
980 : !> \brief clears the buffer
981 : !> \param diis_buffer the buffer to clear
982 : ! **************************************************************************************************
983 3210 : PURE SUBROUTINE qs_diis_b_clear_kp(diis_buffer)
984 :
985 : TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
986 :
987 3210 : diis_buffer%ncall = 0
988 :
989 3210 : END SUBROUTINE qs_diis_b_clear_kp
990 :
991 : ! **************************************************************************************************
992 : !> \brief Update info about the current buffer step ib and the current number of buffers nb
993 : !> \param diis_buffer ...
994 : !> \param ib ...
995 : !> \param nb ...
996 : ! **************************************************************************************************
997 27450 : SUBROUTINE qs_diis_b_info_kp(diis_buffer, ib, nb)
998 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
999 : INTEGER, INTENT(OUT) :: ib, nb
1000 :
1001 27450 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
1002 27450 : diis_buffer%ncall = diis_buffer%ncall + 1
1003 27450 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
1004 :
1005 27450 : END SUBROUTINE qs_diis_b_info_kp
1006 :
1007 : ! **************************************************************************************************
1008 : !> \brief Calculate and store the error for a given k-point
1009 : !> \param diis_buffer ...
1010 : !> \param ib ...
1011 : !> \param mos ...
1012 : !> \param kc ...
1013 : !> \param sc ...
1014 : !> \param ispin ...
1015 : !> \param ikp ...
1016 : !> \param cmos caller-owned scratch; the controller preallocates it and the DIIS buffer
1017 : !> \param real_wfn ...
1018 : !> \param local use local GEMM on a one-rank k-point group
1019 : !> \param gemm_ctx private local GEMM context, unused for distributed groups
1020 : !> \note We assume that we always have an overlap matrix.
1021 : !> TODO: do we need to pass the kp weight for the back Fourier transform?
1022 : ! **************************************************************************************************
1023 97152 : SUBROUTINE qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, cmos, real_wfn, local, gemm_ctx)
1024 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1025 : INTEGER, INTENT(IN) :: ib
1026 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1027 : TYPE(cp_cfm_type), INTENT(INOUT) :: kc, sc
1028 : INTEGER, INTENT(IN) :: ispin, ikp
1029 : TYPE(cp_cfm_type), INTENT(INOUT) :: cmos
1030 : LOGICAL, INTENT(IN), OPTIONAL :: real_wfn, local
1031 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
1032 :
1033 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_calc_err_kp'
1034 :
1035 : INTEGER :: handle, homo, nao
1036 : LOGICAL :: my_real_wfn
1037 : REAL(dp) :: maxocc
1038 : TYPE(cp_cfm_type), POINTER :: new_errors, parameters, smat
1039 : TYPE(cp_fm_type), POINTER :: rmos
1040 :
1041 48576 : NULLIFY (rmos, parameters, new_errors, smat)
1042 :
1043 48576 : CALL timeset(routineN, handle)
1044 48576 : my_real_wfn = .FALSE.
1045 48576 : IF (PRESENT(real_wfn)) my_real_wfn = real_wfn
1046 :
1047 : !Calculate the error for this given k-point, store the KS matrix as well as the ovlp matrix
1048 : !All of this happens within the kp subgroups
1049 :
1050 : ! Quick return, if no DIIS is requested
1051 48576 : IF (diis_buffer%nbuffer < 1) THEN
1052 0 : CALL timestop(handle)
1053 0 : RETURN
1054 : END IF
1055 48576 : CALL get_mo_set(mos(ispin), nao=nao, homo=homo, maxocc=maxocc)
1056 :
1057 48576 : IF (my_real_wfn) THEN
1058 0 : CALL get_mo_set(mos(ispin), mo_coeff=rmos)
1059 : CALL cp_cfm_scale_and_add_fm(CMPLX(0.0_dp, KIND=dp), cmos, &
1060 0 : CMPLX(1.0_dp, KIND=dp), rmos)
1061 : ELSE
1062 48576 : CALL cp_cfm_to_cfm(mos(ispin)%cmo_coeff, cmos)
1063 : END IF
1064 :
1065 48576 : new_errors => diis_buffer%error(ib, ispin, ikp)
1066 48576 : parameters => diis_buffer%param(ib, ispin, ikp)
1067 48576 : smat => diis_buffer%smat(ikp)
1068 :
1069 : !copy the KS and overlap matrices to the DIIS buffer
1070 48576 : CALL cp_cfm_to_cfm(kc, parameters)
1071 48576 : CALL cp_cfm_to_cfm(sc, smat)
1072 :
1073 : ! KC <- K*C; SC <- S*C. Retain the existing factor of two in SC.
1074 : ! new_errors <- KC*(SC)^H - SC*(KC)^H.
1075 48576 : IF (local) THEN
1076 100902 : CPASSERT(PRODUCT(parameters%matrix_struct%context%num_pe) == 1)
1077 100902 : CPASSERT(PRODUCT(cmos%matrix_struct%context%num_pe) == 1)
1078 100902 : CPASSERT(PRODUCT(kc%matrix_struct%context%num_pe) == 1)
1079 100902 : CPASSERT(PRODUCT(smat%matrix_struct%context%num_pe) == 1)
1080 100902 : CPASSERT(PRODUCT(sc%matrix_struct%context%num_pe) == 1)
1081 100902 : CPASSERT(PRODUCT(new_errors%matrix_struct%context%num_pe) == 1)
1082 : CALL gemm_ctx%gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), &
1083 : parameters%local_data, SIZE(parameters%local_data, 1), &
1084 : cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
1085 33634 : kc%local_data, SIZE(kc%local_data, 1))
1086 : CALL gemm_ctx%gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), &
1087 : smat%local_data, SIZE(smat%local_data, 1), &
1088 : cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
1089 33634 : sc%local_data, SIZE(sc%local_data, 1))
1090 : CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
1091 : sc%local_data, SIZE(sc%local_data, 1), &
1092 : kc%local_data, SIZE(kc%local_data, 1), (0.0_dp, 0.0_dp), &
1093 33634 : new_errors%local_data, SIZE(new_errors%local_data, 1))
1094 : CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
1095 : kc%local_data, SIZE(kc%local_data, 1), &
1096 : sc%local_data, SIZE(sc%local_data, 1), (-1.0_dp, 0.0_dp), &
1097 33634 : new_errors%local_data, SIZE(new_errors%local_data, 1))
1098 : ELSE
1099 14942 : CALL parallel_gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), parameters, cmos, (0.0_dp, 0.0_dp), kc)
1100 14942 : CALL parallel_gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), smat, cmos, (0.0_dp, 0.0_dp), sc)
1101 14942 : CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), sc, kc, (0.0_dp, 0.0_dp), new_errors)
1102 14942 : CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), kc, sc, (-1.0_dp, 0.0_dp), new_errors)
1103 : END IF
1104 :
1105 48576 : CALL timestop(handle)
1106 :
1107 : END SUBROUTINE qs_diis_b_calc_err_kp
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points
1111 : !> \param diis_buffer ...
1112 : !> \param coeffs ...
1113 : !> \param ib ...
1114 : !> \param nb ...
1115 : !> \param delta ...
1116 : !> \param error_max ...
1117 : !> \param diis_step ...
1118 : !> \param eps_diis ...
1119 : !> \param nspin ...
1120 : !> \param nkp ...
1121 : !> \param nkp_local ...
1122 : !> \param nmixing ...
1123 : !> \param scf_section ...
1124 : !> \param para_env_inter_kp communicator connecting the k-point groups
1125 : ! **************************************************************************************************
1126 27450 : SUBROUTINE qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, &
1127 : nspin, nkp, nkp_local, nmixing, scf_section, para_env_inter_kp)
1128 :
1129 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1130 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeffs
1131 : INTEGER, INTENT(IN) :: ib, nb
1132 : REAL(KIND=dp), INTENT(IN) :: delta
1133 : REAL(KIND=dp), INTENT(OUT) :: error_max
1134 : LOGICAL, INTENT(OUT) :: diis_step
1135 : REAL(KIND=dp), INTENT(IN) :: eps_diis
1136 : INTEGER, INTENT(IN) :: nspin, nkp, nkp_local
1137 : INTEGER, INTENT(IN), OPTIONAL :: nmixing
1138 : TYPE(section_vals_type), POINTER :: scf_section
1139 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1140 :
1141 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_kp'
1142 : REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
1143 :
1144 : CHARACTER(LEN=2*default_string_length) :: message
1145 : COMPLEX(KIND=dp) :: tmp
1146 27450 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
1147 : INTEGER :: handle, ikp, ispin, jb, my_nmixing, nb1, &
1148 : output_unit
1149 : LOGICAL :: eigenvectors_discarded
1150 27450 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
1151 : TYPE(cp_cfm_type), POINTER :: new_errors
1152 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1153 : TYPE(cp_logger_type), POINTER :: logger
1154 :
1155 27450 : NULLIFY (matrix_struct, new_errors, logger)
1156 :
1157 27450 : CALL timeset(routineN, handle)
1158 :
1159 27450 : diis_step = .FALSE.
1160 :
1161 27450 : my_nmixing = 2
1162 27450 : IF (PRESENT(nmixing)) my_nmixing = nmixing
1163 :
1164 27450 : logger => cp_get_default_logger()
1165 :
1166 : ! Quick return, if no DIIS is requested
1167 27450 : IF (diis_buffer%nbuffer < 1) THEN
1168 0 : CALL timestop(handle)
1169 0 : RETURN
1170 : END IF
1171 :
1172 : ! Check, if a DIIS step is appropriate
1173 27450 : diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
1174 :
1175 : ! Calculate the DIIS buffer, and update it if max_error < eps_diis
1176 27450 : CALL cp_cfm_get_info(diis_buffer%error(ib, 1, 1), matrix_struct=matrix_struct)
1177 109800 : ALLOCATE (b(nb, nb))
1178 27450 : b = 0.0_dp
1179 122578 : DO jb = 1, nb
1180 269808 : DO ikp = 1, nkp_local
1181 405349 : DO ispin = 1, nspin
1182 162991 : new_errors => diis_buffer%error(ib, ispin, ikp)
1183 : CALL cp_cfm_trace(diis_buffer%error(jb, ispin, ikp), new_errors, tmp, &
1184 162991 : conjugate_a=.TRUE., local=.TRUE.)
1185 310221 : b(jb, ib) = b(jb, ib) + 1.0_dp/REAL(nkp, dp)*tmp
1186 : END DO
1187 : END DO
1188 : END DO
1189 : ! Accumulate local AO entries and k points first, then reduce the new history column.
1190 : ! Each inter-group communicator contains one corresponding rank per k-point group.
1191 27450 : CALL matrix_struct%para_env%sum(b(:, ib))
1192 27450 : CALL para_env_inter_kp%sum(b(:, ib))
1193 122578 : DO jb = 1, nb
1194 122578 : b(ib, jb) = CONJG(b(jb, ib))
1195 : END DO
1196 :
1197 27450 : error_max = SQRT(REAL(b(ib, ib))**2 + AIMAG(b(ib, ib))**2)
1198 :
1199 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
1200 27450 : extension=".scfLog")
1201 27450 : IF (output_unit > 0) THEN
1202 : WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
1203 0 : "DIIS | Current SCF DIIS buffer size: ", nb, &
1204 0 : "DIIS | Maximum SCF DIIS error at last step: ", error_max, &
1205 0 : "DIIS | Current SCF convergence: ", delta, &
1206 0 : "DIIS | Threshold value for a DIIS step: ", eps_diis
1207 0 : IF (error_max < eps_diis) THEN
1208 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1209 0 : "DIIS | => The SCF DIIS buffer will be updated"
1210 : ELSE
1211 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1212 0 : "DIIS | => No update of the SCF DIIS buffer"
1213 : END IF
1214 0 : IF (diis_step .AND. (error_max < eps_diis)) THEN
1215 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1216 0 : "DIIS | => A SCF DIIS step will be performed"
1217 : ELSE
1218 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1219 0 : "DIIS | => No SCF DIIS step will be performed"
1220 : END IF
1221 : END IF
1222 :
1223 : ! Update the SCF DIIS buffer
1224 27450 : IF (error_max < eps_diis) THEN
1225 82070 : DO jb = 1, nb
1226 64084 : diis_buffer%b_matrix(ib, jb) = b(ib, jb)
1227 82070 : diis_buffer%b_matrix(jb, ib) = b(jb, ib)
1228 : END DO
1229 : ELSE
1230 :
1231 9464 : diis_step = .FALSE.
1232 : END IF
1233 27450 : DEALLOCATE (b)
1234 :
1235 : ! Perform DIIS step
1236 27450 : IF (diis_step) THEN
1237 :
1238 10652 : nb1 = nb + 1
1239 :
1240 42608 : ALLOCATE (a(nb1, nb1))
1241 31956 : ALLOCATE (b(nb1, nb1))
1242 31956 : ALLOCATE (ev(nb1))
1243 :
1244 : ! Set up the linear DIIS equation system
1245 212932 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
1246 :
1247 51794 : b(1:nb, nb1) = -1.0_dp
1248 51794 : b(nb1, 1:nb) = -1.0_dp
1249 10652 : b(nb1, nb1) = 0.0_dp
1250 :
1251 : ! Solve the linear DIIS equation system
1252 62446 : ev(1:nb1) = 0.0_dp !eigenvalues
1253 316520 : a(1:nb1, 1:nb1) = 0.0_dp !eigenvectors
1254 10652 : CALL diag_complex(b(1:nb1, 1:nb1), a(1:nb1, 1:nb1), ev(1:nb1))
1255 316520 : b(1:nb1, 1:nb1) = a(1:nb1, 1:nb1)
1256 :
1257 10652 : eigenvectors_discarded = .FALSE.
1258 :
1259 62446 : DO jb = 1, nb1
1260 62446 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
1261 18158 : IF (output_unit > 0) THEN
1262 0 : IF (.NOT. eigenvectors_discarded) THEN
1263 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1264 0 : "DIIS | Checking eigenvalues of the DIIS error matrix"
1265 : END IF
1266 : WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
1267 0 : "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
1268 0 : "threshold ", eigenvalue_threshold
1269 0 : CALL compress(message)
1270 0 : WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
1271 0 : eigenvectors_discarded = .TRUE.
1272 : END IF
1273 108664 : a(1:nb1, jb) = 0.0_dp
1274 : ELSE
1275 197204 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
1276 : END IF
1277 : END DO
1278 :
1279 10652 : IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
1280 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1281 0 : "DIIS | The corresponding eigenvectors were discarded"
1282 : END IF
1283 :
1284 327172 : coeffs(1:nb) = -MATMUL(a(1:nb, 1:nb1), CONJG(b(nb1, 1:nb1)))
1285 : ELSE
1286 :
1287 70784 : coeffs(:) = 0.0_dp
1288 16798 : coeffs(ib) = 1.0_dp
1289 : END IF
1290 :
1291 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1292 27450 : "PRINT%DIIS_INFO")
1293 :
1294 27450 : CALL timestop(handle)
1295 :
1296 54900 : END SUBROUTINE qs_diis_b_step_kp
1297 10652 : END MODULE qs_diis
|
