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