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