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