(git:ccc2433)
qs_diis.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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_operations, ONLY: copy_dbcsr_to_fm
31  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
32  cp_fm_scale,&
33  cp_fm_scale_and_add,&
34  cp_fm_symm,&
35  cp_fm_trace
36  USE cp_fm_struct, ONLY: cp_fm_struct_type
37  USE cp_fm_types, ONLY: cp_fm_create,&
38  cp_fm_get_info,&
39  cp_fm_maxabsval,&
40  cp_fm_release,&
41  cp_fm_set_all,&
42  cp_fm_to_fm,&
43  cp_fm_type
44  USE cp_log_handling, ONLY: cp_get_default_logger,&
45  cp_logger_type,&
46  cp_to_string
47  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
48  cp_print_key_unit_nr
49  USE dbcsr_api, ONLY: &
50  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_maxabs, dbcsr_multiply, &
51  dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_transposed, dbcsr_type
52  USE dm_ls_scf_types, ONLY: ls_scf_env_type
53  USE input_section_types, ONLY: section_vals_type
54  USE kinds, ONLY: default_string_length,&
55  dp
56  USE mathlib, ONLY: complex_diag,&
57  diamat_all
58  USE message_passing, ONLY: mp_para_env_type
59  USE parallel_gemm_api, ONLY: parallel_gemm
60  USE qs_diis_types, ONLY: qs_diis_buffer_type,&
61  qs_diis_buffer_type_kp,&
62  qs_diis_buffer_type_sparse
63  USE qs_environment_types, ONLY: get_qs_env,&
64  qs_environment_type
65  USE qs_mo_types, ONLY: get_mo_set,&
66  mo_set_type
67  USE string_utilities, ONLY: compress
68 #include "./base/base_uses.f90"
69 
70  IMPLICIT NONE
71 
72  PRIVATE
73 
74  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_diis'
75 
76  ! Public subroutines
77 
78  PUBLIC :: qs_diis_b_clear, &
79  qs_diis_b_create, &
80  qs_diis_b_step
81  PUBLIC :: qs_diis_b_clear_sparse, &
82  qs_diis_b_create_sparse, &
83  qs_diis_b_step_4lscf
84  PUBLIC :: qs_diis_b_clear_kp, &
85  qs_diis_b_create_kp, &
86  qs_diis_b_step_kp, &
87  qs_diis_b_calc_err_kp, &
88  qs_diis_b_info_kp
89 
90 CONTAINS
91 
92 ! **************************************************************************************************
93 !> \brief Allocates an SCF DIIS buffer
94 !> \param diis_buffer the buffer to create
95 !> \param nbuffer ...
96 !> \par History
97 !> 02.2003 created [fawzi]
98 !> \author fawzi
99 ! **************************************************************************************************
100  SUBROUTINE qs_diis_b_create(diis_buffer, nbuffer)
101 
102  TYPE(qs_diis_buffer_type), INTENT(OUT) :: diis_buffer
103  INTEGER, INTENT(in) :: nbuffer
104 
105  CHARACTER(len=*), PARAMETER :: routinen = 'qs_diis_b_create'
106 
107  INTEGER :: handle
108 
109 ! -------------------------------------------------------------------------
110 
111  CALL timeset(routinen, handle)
112 
113  NULLIFY (diis_buffer%b_matrix)
114  NULLIFY (diis_buffer%error)
115  NULLIFY (diis_buffer%param)
116  diis_buffer%nbuffer = nbuffer
117  diis_buffer%ncall = 0
118 
119  CALL timestop(handle)
120 
121  END SUBROUTINE qs_diis_b_create
122 
123 ! **************************************************************************************************
124 !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
125 !> variables and with a buffer size of nbuffer.
126 !> \param diis_buffer the buffer to initialize
127 !> \param matrix_struct the structure for the matrix of the buffer
128 !> \param nspin ...
129 !> \param scf_section ...
130 !> \par History
131 !> - Creation (07.05.2001, Matthias Krack)
132 !> - Changed to BLACS matrix usage (08.06.2001,MK)
133 !> - DIIS for ROKS (05.04.06,MK)
134 !> \author Matthias Krack
135 !> \note
136 !> check to allocate matrixes only when needed, using a linked list?
137 ! **************************************************************************************************
138  SUBROUTINE qs_diis_b_check_i_alloc(diis_buffer, matrix_struct, nspin, &
139  scf_section)
140 
141  TYPE(qs_diis_buffer_type), INTENT(INOUT) :: diis_buffer
142  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
143  INTEGER, INTENT(IN) :: nspin
144  TYPE(section_vals_type), POINTER :: scf_section
145 
146  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_check_i_alloc'
147 
148  INTEGER :: handle, ibuffer, ispin, nbuffer, &
149  output_unit
150  TYPE(cp_logger_type), POINTER :: logger
151 
152 ! -------------------------------------------------------------------------
153 
154  CALL timeset(routinen, handle)
155 
156  logger => cp_get_default_logger()
157 
158  nbuffer = diis_buffer%nbuffer
159 
160  IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
161  ALLOCATE (diis_buffer%error(nbuffer, nspin))
162 
163  DO ispin = 1, nspin
164  DO ibuffer = 1, nbuffer
165  CALL cp_fm_create(diis_buffer%error(ibuffer, ispin), &
166  name="qs_diis_b%error("// &
167  trim(adjustl(cp_to_string(ibuffer)))//","// &
168  trim(adjustl(cp_to_string(ibuffer)))//")", &
169  matrix_struct=matrix_struct)
170  END DO
171  END DO
172  END IF
173 
174  IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
175  ALLOCATE (diis_buffer%param(nbuffer, nspin))
176 
177  DO ispin = 1, nspin
178  DO ibuffer = 1, nbuffer
179  CALL cp_fm_create(diis_buffer%param(ibuffer, ispin), &
180  name="qs_diis_b%param("// &
181  trim(adjustl(cp_to_string(ibuffer)))//","// &
182  trim(adjustl(cp_to_string(ibuffer)))//")", &
183  matrix_struct=matrix_struct)
184  END DO
185  END DO
186  END IF
187 
188  IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
189  ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
190  diis_buffer%b_matrix = 0.0_dp
191  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
192  extension=".scfLog")
193  IF (output_unit > 0) THEN
194  WRITE (unit=output_unit, fmt="(/,T9,A)") &
195  "DIIS | The SCF DIIS buffer was allocated and initialized"
196  END IF
197  CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
198  "PRINT%DIIS_INFO")
199  END IF
200 
201  CALL timestop(handle)
202 
203  END SUBROUTINE qs_diis_b_check_i_alloc
204 
205 ! **************************************************************************************************
206 !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step.
207 !> \param diis_buffer ...
208 !> \param mo_array ...
209 !> \param kc ...
210 !> \param sc ...
211 !> \param delta ...
212 !> \param error_max ...
213 !> \param diis_step ...
214 !> \param eps_diis ...
215 !> \param nmixing ...
216 !> \param s_matrix ...
217 !> \param scf_section ...
218 !> \param roks ...
219 !> \par History
220 !> - Creation (07.05.2001, Matthias Krack)
221 !> - Changed to BLACS matrix usage (08.06.2001, MK)
222 !> - 03.2003 rewamped [fawzi]
223 !> - Adapted for high-spin ROKS (08.04.06,MK)
224 !> \author Matthias Krack
225 ! **************************************************************************************************
226  SUBROUTINE qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, &
227  diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
228 
229  TYPE(qs_diis_buffer_type), POINTER :: diis_buffer
230  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
231  TYPE(cp_fm_type), DIMENSION(:), POINTER :: kc
232  TYPE(cp_fm_type), INTENT(IN) :: sc
233  REAL(kind=dp), INTENT(IN) :: delta
234  REAL(kind=dp), INTENT(OUT) :: error_max
235  LOGICAL, INTENT(OUT) :: diis_step
236  REAL(kind=dp), INTENT(IN) :: eps_diis
237  INTEGER, INTENT(IN), OPTIONAL :: nmixing
238  TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
239  POINTER :: s_matrix
240  TYPE(section_vals_type), POINTER :: scf_section
241  LOGICAL, INTENT(IN), OPTIONAL :: roks
242 
243  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_step'
244  REAL(kind=dp), PARAMETER :: eigenvalue_threshold = 1.0e-12_dp
245 
246  CHARACTER(LEN=2*default_string_length) :: message
247  INTEGER :: handle, homo, ib, imo, ispin, jb, &
248  my_nmixing, nao, nb, nb1, nmo, nspin, &
249  output_unit
250  LOGICAL :: eigenvectors_discarded, my_roks
251  REAL(kind=dp) :: maxocc, tmp
252  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev, occ
253  REAL(kind=dp), DIMENSION(:), POINTER :: occa, occb
254  REAL(kind=dp), DIMENSION(:, :), POINTER :: a, b
255  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
256  TYPE(cp_fm_type), POINTER :: c, new_errors, old_errors, parameters
257  TYPE(cp_logger_type), POINTER :: logger
258 
259 ! -------------------------------------------------------------------------
260 
261  CALL timeset(routinen, handle)
262 
263  nspin = SIZE(mo_array)
264  diis_step = .false.
265 
266  IF (PRESENT(roks)) THEN
267  my_roks = .true.
268  nspin = 1
269  ELSE
270  my_roks = .false.
271  END IF
272 
273  my_nmixing = 2
274  IF (PRESENT(nmixing)) my_nmixing = nmixing
275 
276  NULLIFY (c, new_errors, old_errors, parameters, matrix_struct, a, b, occa, occb)
277  logger => cp_get_default_logger()
278 
279  ! Quick return, if no DIIS is requested
280 
281  IF (diis_buffer%nbuffer < 1) THEN
282  CALL timestop(handle)
283  RETURN
284  END IF
285 
286  CALL cp_fm_get_info(kc(1), &
287  matrix_struct=matrix_struct)
288  CALL qs_diis_b_check_i_alloc(diis_buffer, &
289  matrix_struct=matrix_struct, &
290  nspin=nspin, &
291  scf_section=scf_section)
292 
293  error_max = 0.0_dp
294 
295  ib = modulo(diis_buffer%ncall, diis_buffer%nbuffer) + 1
296  diis_buffer%ncall = diis_buffer%ncall + 1
297  nb = min(diis_buffer%ncall, diis_buffer%nbuffer)
298 
299  DO ispin = 1, nspin
300 
301  CALL get_mo_set(mo_set=mo_array(ispin), &
302  nao=nao, &
303  nmo=nmo, &
304  homo=homo, &
305  mo_coeff=c, &
306  occupation_numbers=occa, &
307  maxocc=maxocc)
308 
309  new_errors => diis_buffer%error(ib, ispin)
310  parameters => diis_buffer%param(ib, ispin)
311 
312  ! Copy the Kohn-Sham matrix K to the DIIS buffer
313 
314  CALL cp_fm_to_fm(kc(ispin), parameters)
315 
316  IF (my_roks) THEN
317 
318  ALLOCATE (occ(nmo))
319 
320  CALL get_mo_set(mo_set=mo_array(2), &
321  occupation_numbers=occb)
322 
323  DO imo = 1, nmo
324  occ(imo) = sqrt(occa(imo) + occb(imo))
325  END DO
326 
327  CALL cp_fm_to_fm(c, sc)
328  CALL cp_fm_column_scale(sc, occ(1:homo))
329 
330  ! KC <- K*C
331  CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, parameters, sc, 0.0_dp, kc(ispin))
332 
333  IF (PRESENT(s_matrix)) THEN
334  CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
335  ! SC <- S*C
336  CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, new_errors, c, 0.0_dp, sc)
337  CALL cp_fm_column_scale(sc, occ(1:homo))
338  END IF
339 
340  ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
341  ! or for an orthogonal basis
342  ! new_errors <- KC*C^T - C*(KC)^T = K*P - P*K with S = I
343  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
344  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
345 
346  DEALLOCATE (occ)
347 
348  ELSE
349 
350  ! KC <- K*C
351  CALL cp_fm_symm("L", "U", nao, homo, maxocc, parameters, c, 0.0_dp, kc(ispin))
352 
353  IF (PRESENT(s_matrix)) THEN
354  ! I guess that this copy can be avoided for LSD
355  CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
356  ! sc <- S*C
357  CALL cp_fm_symm("L", "U", nao, homo, 2.0_dp, new_errors, c, 0.0_dp, sc)
358  ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
359  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
360  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
361  ELSE
362  ! new_errors <- KC*(C)^T - C*(KC)^T = K*P - P*K
363  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, c, kc(ispin), 0.0_dp, new_errors)
364  CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), c, -1.0_dp, new_errors)
365  END IF
366 
367  END IF
368 
369  CALL cp_fm_maxabsval(new_errors, tmp)
370  error_max = max(error_max, tmp)
371 
372  END DO
373 
374  ! Check, if a DIIS step is appropriate
375 
376  diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
377 
378  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
379  extension=".scfLog")
380  IF (output_unit > 0) THEN
381  WRITE (unit=output_unit, fmt="(/,T9,A,I4,/,(T9,A,ES12.3))") &
382  "DIIS | Current SCF DIIS buffer size: ", nb, &
383  "DIIS | Maximum SCF DIIS error vector element:", error_max, &
384  "DIIS | Current SCF convergence: ", delta, &
385  "DIIS | Threshold value for a DIIS step: ", eps_diis
386  IF (error_max < eps_diis) THEN
387  WRITE (unit=output_unit, fmt="(T9,A)") &
388  "DIIS | => The SCF DIIS buffer will be updated"
389  ELSE
390  WRITE (unit=output_unit, fmt="(T9,A)") &
391  "DIIS | => No update of the SCF DIIS buffer"
392  END IF
393  IF (diis_step .AND. (error_max < eps_diis)) THEN
394  WRITE (unit=output_unit, fmt="(T9,A,/)") &
395  "DIIS | => A SCF DIIS step will be performed"
396  ELSE
397  WRITE (unit=output_unit, fmt="(T9,A,/)") &
398  "DIIS | => No SCF DIIS step will be performed"
399  END IF
400  END IF
401 
402  ! Update the SCF DIIS buffer
403 
404  IF (error_max < eps_diis) THEN
405 
406  b => diis_buffer%b_matrix
407 
408  DO jb = 1, nb
409  b(jb, ib) = 0.0_dp
410  DO ispin = 1, nspin
411  old_errors => diis_buffer%error(jb, ispin)
412  new_errors => diis_buffer%error(ib, ispin)
413  CALL cp_fm_trace(old_errors, new_errors, tmp)
414  b(jb, ib) = b(jb, ib) + tmp
415  END DO
416  b(ib, jb) = b(jb, ib)
417  END DO
418 
419  ELSE
420 
421  diis_step = .false.
422 
423  END IF
424 
425  ! Perform DIIS step
426 
427  IF (diis_step) THEN
428 
429  nb1 = nb + 1
430 
431  ALLOCATE (a(nb1, nb1))
432  ALLOCATE (b(nb1, nb1))
433  ALLOCATE (ev(nb1))
434 
435  ! Set up the linear DIIS equation system
436 
437  b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
438 
439  b(1:nb, nb1) = -1.0_dp
440  b(nb1, 1:nb) = -1.0_dp
441  b(nb1, nb1) = 0.0_dp
442 
443  ! Solve the linear DIIS equation system
444 
445  ev(1:nb1) = 0.0_dp
446  CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
447 
448  a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
449 
450  eigenvectors_discarded = .false.
451 
452  DO jb = 1, nb1
453  IF (abs(ev(jb)) < eigenvalue_threshold) THEN
454  IF (output_unit > 0) THEN
455  IF (.NOT. eigenvectors_discarded) THEN
456  WRITE (unit=output_unit, fmt="(T9,A)") &
457  "DIIS | Checking eigenvalues of the DIIS error matrix"
458  END IF
459  WRITE (unit=message, fmt="(T9,A,I6,A,ES10.1,A,ES10.1)") &
460  "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
461  "threshold ", eigenvalue_threshold
462  CALL compress(message)
463  WRITE (unit=output_unit, fmt="(T9,A)") trim(message)
464  eigenvectors_discarded = .true.
465  END IF
466  a(1:nb1, jb) = 0.0_dp
467  ELSE
468  a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
469  END IF
470  END DO
471 
472  IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
473  WRITE (unit=output_unit, fmt="(T9,A,/)") &
474  "DIIS | The corresponding eigenvectors were discarded"
475  END IF
476 
477  ev(1:nb) = matmul(a(1:nb, 1:nb1), b(nb1, 1:nb1))
478 
479  ! Update Kohn-Sham matrix
480 
481  DO ispin = 1, nspin
482  CALL cp_fm_set_all(kc(ispin), 0.0_dp)
483  DO jb = 1, nb
484  parameters => diis_buffer%param(jb, ispin)
485  CALL cp_fm_scale_and_add(1.0_dp, kc(ispin), -ev(jb), parameters)
486  END DO
487  END DO
488 
489  DEALLOCATE (a)
490  DEALLOCATE (b)
491  DEALLOCATE (ev)
492 
493  ELSE
494 
495  DO ispin = 1, nspin
496  parameters => diis_buffer%param(ib, ispin)
497  CALL cp_fm_to_fm(parameters, kc(ispin))
498  END DO
499 
500  END IF
501 
502  CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
503  "PRINT%DIIS_INFO")
504 
505  CALL timestop(handle)
506 
507  END SUBROUTINE qs_diis_b_step
508 
509 ! **************************************************************************************************
510 !> \brief clears the buffer
511 !> \param diis_buffer the buffer to clear
512 !> \par History
513 !> 02.2003 created [fawzi]
514 !> \author fawzi
515 ! **************************************************************************************************
516  PURE SUBROUTINE qs_diis_b_clear(diis_buffer)
517 
518  TYPE(qs_diis_buffer_type), INTENT(INOUT) :: diis_buffer
519 
520  diis_buffer%ncall = 0
521 
522  END SUBROUTINE qs_diis_b_clear
523 
524 ! **************************************************************************************************
525 !> \brief Update the SCF DIIS buffer in linear scaling SCF (LS-SCF),
526 !> and if appropriate does a diis step.
527 !> \param diis_buffer ...
528 !> \param qs_env ...
529 !> \param ls_scf_env ...
530 !> \param unit_nr ...
531 !> \param iscf ...
532 !> \param diis_step ...
533 !> \param eps_diis ...
534 !> \param nmixing ...
535 !> \param s_matrix ...
536 !> \param threshold ...
537 !> \par History
538 !> - Adapted for LS-SCF (10-11-14) from qs_diis_b_step
539 !> \author Fredy W. Aquino
540 ! **************************************************************************************************
541 
542  SUBROUTINE qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, &
543  diis_step, eps_diis, nmixing, s_matrix, threshold)
544 ! Note.- Input: ls_scf_env%matrix_p(ispin) , Density Matrix
545 ! matrix_ks (from qs_env) , Kohn-Sham Matrix (IN/OUT)
546 
547  TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer
548  TYPE(qs_environment_type), POINTER :: qs_env
549  TYPE(ls_scf_env_type) :: ls_scf_env
550  INTEGER, INTENT(IN) :: unit_nr, iscf
551  LOGICAL, INTENT(OUT) :: diis_step
552  REAL(kind=dp), INTENT(IN) :: eps_diis
553  INTEGER, INTENT(IN), OPTIONAL :: nmixing
554  TYPE(dbcsr_type), OPTIONAL :: s_matrix
555  REAL(kind=dp), INTENT(IN) :: threshold
556 
557  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_step_4lscf'
558  REAL(kind=dp), PARAMETER :: eigenvalue_threshold = 1.0e-12_dp
559 
560  INTEGER :: handle, ib, ispin, jb, my_nmixing, nb, &
561  nb1, nspin
562  REAL(kind=dp) :: error_max, tmp
563  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev
564  REAL(kind=dp), DIMENSION(:, :), POINTER :: a, b
565  TYPE(cp_logger_type), POINTER :: logger
566  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
567  TYPE(dbcsr_type) :: matrix_kserr_t, matrix_tmp
568  TYPE(dbcsr_type), POINTER :: new_errors, old_errors, parameters
569  TYPE(mp_para_env_type), POINTER :: para_env
570 
571  CALL timeset(routinen, handle)
572  nspin = ls_scf_env%nspins
573  diis_step = .false.
574  my_nmixing = 2
575  IF (PRESENT(nmixing)) my_nmixing = nmixing
576  NULLIFY (new_errors, old_errors, parameters, a, b)
577  logger => cp_get_default_logger()
578  ! Quick return, if no DIIS is requested
579  IF (diis_buffer%nbuffer < 1) THEN
580  CALL timestop(handle)
581  RETURN
582  END IF
583 
584 ! Getting current Kohn-Sham matrix from qs_env
585  CALL get_qs_env(qs_env, &
586  para_env=para_env, &
587  matrix_ks=matrix_ks)
588  CALL qs_diis_b_check_i_alloc_sparse( &
589  diis_buffer, &
590  ls_scf_env, &
591  nspin)
592  error_max = 0.0_dp
593 
594  ib = modulo(diis_buffer%ncall, diis_buffer%nbuffer) + 1
595  diis_buffer%ncall = diis_buffer%ncall + 1
596  nb = min(diis_buffer%ncall, diis_buffer%nbuffer)
597 ! Create scratch arrays
598  CALL dbcsr_create(matrix_tmp, &
599  template=ls_scf_env%matrix_ks(1), &
600  matrix_type='N')
601  CALL dbcsr_set(matrix_tmp, 0.0_dp) ! reset matrix
602  CALL dbcsr_create(matrix_kserr_t, &
603  template=ls_scf_env%matrix_ks(1), &
604  matrix_type='N')
605  CALL dbcsr_set(matrix_kserr_t, 0.0_dp) ! reset matrix
606 
607  DO ispin = 1, nspin ! ------ Loop-ispin----START
608 
609  new_errors => diis_buffer%error(ib, ispin)%matrix
610  parameters => diis_buffer%param(ib, ispin)%matrix
611  ! Copy the Kohn-Sham matrix K to the DIIS buffer
612  CALL dbcsr_copy(parameters, & ! out
613  matrix_ks(ispin)%matrix) ! in
614 
615  IF (PRESENT(s_matrix)) THEN ! if-s_matrix ---------- START
616 ! Calculate Kohn-Sham error (non-orthogonal)= K*P*S-(K*P*S)^T
617 ! matrix_tmp = P*S
618  CALL dbcsr_multiply("N", "N", &
619  1.0_dp, ls_scf_env%matrix_p(ispin), &
620  s_matrix, &
621  0.0_dp, matrix_tmp, &
622  filter_eps=threshold)
623 ! new_errors= K*P*S
624  CALL dbcsr_multiply("N", "N", &
625  1.0_dp, matrix_ks(ispin)%matrix, &
626  matrix_tmp, &
627  0.0_dp, new_errors, &
628  filter_eps=threshold)
629 ! matrix_KSerr_t= transpose(K*P*S)
630  CALL dbcsr_transposed(matrix_kserr_t, &
631  new_errors)
632 ! new_errors=K*P*S-transpose(K*P*S)
633  CALL dbcsr_add(new_errors, &
634  matrix_kserr_t, &
635  1.0_dp, -1.0_dp)
636  ELSE ! if-s_matrix ---------- MID
637 ! Calculate Kohn-Sham error (orthogonal)= K*P - P*K
638 ! new_errors=K*P
639  CALL dbcsr_multiply("N", "N", &
640  1.0_dp, matrix_ks(ispin)%matrix, &
641  ls_scf_env%matrix_p(ispin), &
642  0.0_dp, new_errors, &
643  filter_eps=threshold)
644 ! matrix_KSerr_t= transpose(K*P)
645  CALL dbcsr_transposed(matrix_kserr_t, &
646  new_errors)
647 ! new_errors=K*P-transpose(K*P)
648  CALL dbcsr_add(new_errors, &
649  matrix_kserr_t, &
650  1.0_dp, -1.0_dp)
651  END IF ! if-s_matrix ---------- END
652 
653  tmp = dbcsr_maxabs(new_errors)
654  error_max = max(error_max, tmp)
655 
656  END DO ! ------ Loop-ispin----END
657 
658  ! Check, if a DIIS step is appropriate
659 
660  diis_step = (diis_buffer%ncall >= my_nmixing)
661 
662  IF (unit_nr > 0) THEN
663  WRITE (unit_nr, '(A29,I3,A3,4(I3,A1))') &
664  "DIIS: (ncall,nbuffer,ib,nb)=(", iscf, ")=(", &
665  diis_buffer%ncall, ",", diis_buffer%nbuffer, ",", ib, ",", nb, ")"
666  WRITE (unit_nr, '(A57,I3,A3,L1,A1,F10.8,A1,F4.2,A1,L1,A1)') &
667  "DIIS: (diis_step,error_max,eps_diis,error_max<eps_diis)=(", &
668  iscf, ")=(", diis_step, ",", error_max, ",", eps_diis, ",", &
669  (error_max < eps_diis), ")"
670  WRITE (unit_nr, '(A75)') &
671  "DIIS: diis_step=T : Perform DIIS error_max<eps_diis=T : Update DIIS buffer"
672  END IF
673 
674  ! Update the SCF DIIS buffer
675  IF (error_max < eps_diis) THEN
676  b => diis_buffer%b_matrix
677  DO jb = 1, nb
678  b(jb, ib) = 0.0_dp
679  DO ispin = 1, nspin
680  old_errors => diis_buffer%error(jb, ispin)%matrix
681  new_errors => diis_buffer%error(ib, ispin)%matrix
682  CALL dbcsr_dot(old_errors, &
683  new_errors, &
684  tmp) ! out : < f_i | f_j >
685  b(jb, ib) = b(jb, ib) + tmp
686  END DO ! end-loop-ispin
687  b(ib, jb) = b(jb, ib)
688  END DO ! end-loop-jb
689  ELSE
690  diis_step = .false.
691  END IF
692 
693  ! Perform DIIS step
694  IF (diis_step) THEN
695  nb1 = nb + 1
696  ALLOCATE (a(nb1, nb1))
697  ALLOCATE (b(nb1, nb1))
698  ALLOCATE (ev(nb1))
699  ! Set up the linear DIIS equation system
700  b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
701  b(1:nb, nb1) = -1.0_dp
702  b(nb1, 1:nb) = -1.0_dp
703  b(nb1, nb1) = 0.0_dp
704  ! Solve the linear DIIS equation system
705  CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
706  a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
707  DO jb = 1, nb1
708  IF (abs(ev(jb)) < eigenvalue_threshold) THEN
709  a(1:nb1, jb) = 0.0_dp
710  ELSE
711  a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
712  END IF
713  END DO ! end-loop-jb
714 
715  ev(1:nb) = matmul(a(1:nb, 1:nb1), b(nb1, 1:nb1))
716 
717  ! Update Kohn-Sham matrix
718  IF (iscf .GE. ls_scf_env%iter_ini_diis) THEN ! if-iscf-to-updateKS------ START
719 
720  IF (unit_nr > 0) THEN
721  WRITE (unit_nr, '(A40,I3)') 'DIIS: Updating Kohn-Sham matrix at iscf=', iscf
722  END IF
723 
724  DO ispin = 1, nspin
725  CALL dbcsr_set(matrix_ks(ispin)%matrix, & ! reset matrix
726  0.0_dp)
727  DO jb = 1, nb
728  parameters => diis_buffer%param(jb, ispin)%matrix
729  CALL dbcsr_add(matrix_ks(ispin)%matrix, parameters, &
730  1.0_dp, -ev(jb))
731  END DO ! end-loop-jb
732  END DO ! end-loop-ispin
733  END IF ! if-iscf-to-updateKS------ END
734 
735  DEALLOCATE (a)
736  DEALLOCATE (b)
737  DEALLOCATE (ev)
738 
739  ELSE
740  DO ispin = 1, nspin
741  parameters => diis_buffer%param(ib, ispin)%matrix
742  CALL dbcsr_copy(parameters, & ! out
743  matrix_ks(ispin)%matrix) ! in
744  END DO ! end-loop-ispin
745  END IF
746  CALL dbcsr_release(matrix_tmp)
747  CALL dbcsr_release(matrix_kserr_t)
748  CALL timestop(handle)
749 
750  END SUBROUTINE qs_diis_b_step_4lscf
751 
752 ! **************************************************************************************************
753 !> \brief Allocate and initialize a DIIS buffer with a buffer size of nbuffer.
754 !> \param diis_buffer the buffer to initialize
755 !> \param ls_scf_env ...
756 !> \param nspin ...
757 !> \par History
758 !> - Adapted from qs_diis_b_check_i_alloc for sparse matrices and
759 !> used in LS-SCF module (ls_scf_main) (10-11-14)
760 !> \author Fredy W. Aquino
761 !> \note
762 !> check to allocate matrices only when needed
763 ! **************************************************************************************************
764 
765  SUBROUTINE qs_diis_b_check_i_alloc_sparse(diis_buffer, ls_scf_env, &
766  nspin)
767 
768  TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
769  TYPE(ls_scf_env_type) :: ls_scf_env
770  INTEGER, INTENT(IN) :: nspin
771 
772  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_check_i_alloc_sparse'
773 
774  INTEGER :: handle, ibuffer, ispin, nbuffer
775  TYPE(cp_logger_type), POINTER :: logger
776 
777 ! -------------------------------------------------------------------------
778 
779  CALL timeset(routinen, handle)
780 
781  logger => cp_get_default_logger()
782 
783  nbuffer = diis_buffer%nbuffer
784 
785  IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
786  ALLOCATE (diis_buffer%error(nbuffer, nspin))
787 
788  DO ispin = 1, nspin
789  DO ibuffer = 1, nbuffer
790  ALLOCATE (diis_buffer%error(ibuffer, ispin)%matrix)
791 
792  CALL dbcsr_create(diis_buffer%error(ibuffer, ispin)%matrix, &
793  template=ls_scf_env%matrix_ks(1), &
794  matrix_type='N')
795  END DO
796  END DO
797  END IF
798 
799  IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
800  ALLOCATE (diis_buffer%param(nbuffer, nspin))
801 
802  DO ispin = 1, nspin
803  DO ibuffer = 1, nbuffer
804  ALLOCATE (diis_buffer%param(ibuffer, ispin)%matrix)
805  CALL dbcsr_create(diis_buffer%param(ibuffer, ispin)%matrix, &
806  template=ls_scf_env%matrix_ks(1), &
807  matrix_type='N')
808  END DO
809  END DO
810  END IF
811 
812  IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
813  ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
814 
815  diis_buffer%b_matrix = 0.0_dp
816  END IF
817 
818  CALL timestop(handle)
819 
820  END SUBROUTINE qs_diis_b_check_i_alloc_sparse
821 
822 ! **************************************************************************************************
823 !> \brief clears the DIIS buffer in LS-SCF calculation
824 !> \param diis_buffer the buffer to clear
825 !> \par History
826 !> 10-11-14 created [FA] modified from qs_diis_b_clear
827 !> \author Fredy W. Aquino
828 ! **************************************************************************************************
829 
830  PURE SUBROUTINE qs_diis_b_clear_sparse(diis_buffer)
831 
832  TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
833 
834  diis_buffer%ncall = 0
835 
836  END SUBROUTINE qs_diis_b_clear_sparse
837 
838 ! **************************************************************************************************
839 !> \brief Allocates an SCF DIIS buffer for LS-SCF calculation
840 !> \param diis_buffer the buffer to create
841 !> \param nbuffer ...
842 !> \par History
843 !> 10-11-14 created [FA] modified from qs_diis_b_create
844 !> \author Fredy W. Aquino
845 ! **************************************************************************************************
846  PURE SUBROUTINE qs_diis_b_create_sparse(diis_buffer, nbuffer)
847 
848  TYPE(qs_diis_buffer_type_sparse), INTENT(OUT) :: diis_buffer
849  INTEGER, INTENT(in) :: nbuffer
850 
851  NULLIFY (diis_buffer%b_matrix)
852  NULLIFY (diis_buffer%error)
853  NULLIFY (diis_buffer%param)
854  diis_buffer%nbuffer = nbuffer
855  diis_buffer%ncall = 0
856 
857  END SUBROUTINE qs_diis_b_create_sparse
858 
859 ! **************************************************************************************************
860 !> \brief Allocates an SCF DIIS buffer for k-points
861 !> \param diis_buffer the buffer to create
862 !> \param nbuffer ...
863 ! **************************************************************************************************
864  SUBROUTINE qs_diis_b_create_kp(diis_buffer, nbuffer)
865 
866  TYPE(qs_diis_buffer_type_kp), INTENT(OUT) :: diis_buffer
867  INTEGER, INTENT(in) :: nbuffer
868 
869  NULLIFY (diis_buffer%b_matrix)
870  NULLIFY (diis_buffer%error)
871  NULLIFY (diis_buffer%param)
872  NULLIFY (diis_buffer%smat)
873  diis_buffer%nbuffer = nbuffer
874  diis_buffer%ncall = 0
875 
876  END SUBROUTINE qs_diis_b_create_kp
877 
878 ! **************************************************************************************************
879 !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
880 !> variables and with a buffer size of nbuffer, in the k-point case
881 !> \param diis_buffer the buffer to initialize
882 !> \param matrix_struct the structure for the matrix of the buffer note: this is in the kp subgroup
883 !> \param nspin ...
884 !> \param nkp ...
885 !> \param scf_section ...
886 ! **************************************************************************************************
887  SUBROUTINE qs_diis_b_check_i_alloc_kp(diis_buffer, matrix_struct, nspin, nkp, scf_section)
888 
889  TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
890  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
891  INTEGER, INTENT(IN) :: nspin, nkp
892  TYPE(section_vals_type), POINTER :: scf_section
893 
894  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_check_i_alloc_kp'
895 
896  INTEGER :: handle, ibuffer, ikp, ispin, nbuffer, &
897  output_unit
898  TYPE(cp_logger_type), POINTER :: logger
899 
900 ! -------------------------------------------------------------------------
901 
902  CALL timeset(routinen, handle)
903 
904  logger => cp_get_default_logger()
905 
906  nbuffer = diis_buffer%nbuffer
907 
908  IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
909  ALLOCATE (diis_buffer%error(nbuffer, nspin, nkp))
910 
911  DO ikp = 1, nkp
912  DO ispin = 1, nspin
913  DO ibuffer = 1, nbuffer
914  CALL cp_cfm_create(diis_buffer%error(ibuffer, ispin, ikp), &
915  name="qs_diis_b%error("// &
916  trim(adjustl(cp_to_string(ibuffer)))//","// &
917  trim(adjustl(cp_to_string(ibuffer)))//")", &
918  matrix_struct=matrix_struct)
919  END DO
920  END DO
921  END DO
922  END IF
923 
924  IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
925  ALLOCATE (diis_buffer%param(nbuffer, nspin, nkp))
926 
927  DO ikp = 1, nkp
928  DO ispin = 1, nspin
929  DO ibuffer = 1, nbuffer
930  CALL cp_cfm_create(diis_buffer%param(ibuffer, ispin, ikp), &
931  name="qs_diis_b%param("// &
932  trim(adjustl(cp_to_string(ibuffer)))//","// &
933  trim(adjustl(cp_to_string(ibuffer)))//")", &
934  matrix_struct=matrix_struct)
935  END DO
936  END DO
937  END DO
938  END IF
939 
940  IF (.NOT. ASSOCIATED(diis_buffer%smat)) THEN
941  ALLOCATE (diis_buffer%smat(nkp))
942  DO ikp = 1, nkp
943  CALL cp_cfm_create(diis_buffer%smat(ikp), &
944  name="kp_cfm_smat("// &
945  trim(adjustl(cp_to_string(ibuffer)))//","// &
946  trim(adjustl(cp_to_string(ibuffer)))//")", &
947  matrix_struct=matrix_struct)
948  END DO
949  END IF
950 
951  IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
952  ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
953  diis_buffer%b_matrix = 0.0_dp
954  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
955  extension=".scfLog")
956  IF (output_unit > 0) THEN
957  WRITE (unit=output_unit, fmt="(/,T9,A)") &
958  "DIIS | The SCF DIIS buffer was allocated and initialized"
959  END IF
960  CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
961  "PRINT%DIIS_INFO")
962  END IF
963 
964  CALL timestop(handle)
965 
966  END SUBROUTINE qs_diis_b_check_i_alloc_kp
967 
968 ! **************************************************************************************************
969 !> \brief clears the buffer
970 !> \param diis_buffer the buffer to clear
971 ! **************************************************************************************************
972  PURE SUBROUTINE qs_diis_b_clear_kp(diis_buffer)
973 
974  TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
975 
976  diis_buffer%ncall = 0
977 
978  END SUBROUTINE qs_diis_b_clear_kp
979 
980 ! **************************************************************************************************
981 !> \brief Update info about the current buffer step ib and the current number of buffers nb
982 !> \param diis_buffer ...
983 !> \param ib ...
984 !> \param nb ...
985 ! **************************************************************************************************
986  SUBROUTINE qs_diis_b_info_kp(diis_buffer, ib, nb)
987  TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
988  INTEGER, INTENT(OUT) :: ib, nb
989 
990  ib = modulo(diis_buffer%ncall, diis_buffer%nbuffer) + 1
991  diis_buffer%ncall = diis_buffer%ncall + 1
992  nb = min(diis_buffer%ncall, diis_buffer%nbuffer)
993 
994  END SUBROUTINE qs_diis_b_info_kp
995 
996 ! **************************************************************************************************
997 !> \brief Calculate and store the error for a given k-point
998 !> \param diis_buffer ...
999 !> \param ib ...
1000 !> \param mos ...
1001 !> \param kc ...
1002 !> \param sc ...
1003 !> \param ispin ...
1004 !> \param ikp ...
1005 !> \param nkp_local ...
1006 !> \param scf_section ...
1007 !> \note We assume that we always have an overlap matrix and complex matrices
1008 !> TODO: do we need to pass the kp weight for the back Fourier transform?
1009 ! **************************************************************************************************
1010  SUBROUTINE qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, nkp_local, scf_section)
1011  TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1012  INTEGER, INTENT(IN) :: ib
1013  TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
1014  TYPE(cp_cfm_type), INTENT(INOUT) :: kc, sc
1015  INTEGER, INTENT(IN) :: ispin, ikp, nkp_local
1016  TYPE(section_vals_type), POINTER :: scf_section
1017 
1018  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_calc_err_kp'
1019 
1020  INTEGER :: handle, homo, nao, nmo, nspin
1021  REAL(dp) :: maxocc
1022  TYPE(cp_cfm_type) :: cmos
1023  TYPE(cp_cfm_type), POINTER :: new_errors, parameters, smat
1024  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1025  TYPE(cp_fm_type), POINTER :: imos, rmos
1026 
1027  NULLIFY (matrix_struct, imos, rmos, parameters, new_errors, smat)
1028 
1029  CALL timeset(routinen, handle)
1030 
1031  !Calculate the error for this given k-point, store the KS matrix as well as the ovlp matrix
1032  !All of this happens within the kp subgroups
1033 
1034  ! Quick return, if no DIIS is requested
1035  IF (diis_buffer%nbuffer < 1) THEN
1036  CALL timestop(handle)
1037  RETURN
1038  END IF
1039  nspin = SIZE(mos, 2)
1040 
1041  CALL cp_cfm_get_info(kc, matrix_struct=matrix_struct)
1042  CALL qs_diis_b_check_i_alloc_kp(diis_buffer, &
1043  matrix_struct=matrix_struct, &
1044  nspin=nspin, nkp=nkp_local, &
1045  scf_section=scf_section)
1046 
1047  !We calculate: e(ikp) = F(ikp)*P(ikp)*S(ikp) - S(ikp)*P(ikp)*F(ikp)
1048  CALL get_mo_set(mos(1, ispin), nao=nao, nmo=nmo, homo=homo, mo_coeff=rmos, maxocc=maxocc)
1049  CALL get_mo_set(mos(2, ispin), mo_coeff=imos)
1050  NULLIFY (matrix_struct)
1051  CALL cp_fm_get_info(rmos, matrix_struct=matrix_struct)
1052  CALL cp_cfm_create(cmos, matrix_struct)
1053  CALL cp_fm_to_cfm(rmos, imos, cmos)
1054 
1055  new_errors => diis_buffer%error(ib, ispin, ikp)
1056  parameters => diis_buffer%param(ib, ispin, ikp)
1057  smat => diis_buffer%smat(ikp)
1058 
1059  !copy the KS and overlap matrices to the DIIS buffer
1060  CALL cp_cfm_to_cfm(kc, parameters)
1061  CALL cp_cfm_to_cfm(sc, smat)
1062 
1063  ! KC <- K*C
1064  CALL parallel_gemm("N", "N", nao, homo, nao, cmplx(maxocc, kind=dp), parameters, cmos, (0.0_dp, 0.0_dp), kc)
1065  ! SC <- S*C
1066  CALL parallel_gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), smat, cmos, (0.0_dp, 0.0_dp), sc)
1067 
1068  ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
1069  CALL parallel_gemm("N", "T", nao, nao, homo, (1.0_dp, 0.0_dp), sc, kc, (0.0_dp, 0.0_dp), new_errors)
1070  CALL parallel_gemm("N", "T", nao, nao, homo, (1.0_dp, 0.0_dp), kc, sc, (-1.0_dp, 0.0_dp), new_errors)
1071 
1072  !clean-up
1073  CALL cp_cfm_release(cmos)
1074 
1075  CALL timestop(handle)
1076 
1077  END SUBROUTINE qs_diis_b_calc_err_kp
1078 
1079 ! **************************************************************************************************
1080 !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points
1081 !> \param diis_buffer ...
1082 !> \param coeffs ...
1083 !> \param ib ...
1084 !> \param nb ...
1085 !> \param delta ...
1086 !> \param error_max ...
1087 !> \param diis_step ...
1088 !> \param eps_diis ...
1089 !> \param nspin ...
1090 !> \param nkp ...
1091 !> \param nkp_local ...
1092 !> \param nmixing ...
1093 !> \param scf_section ...
1094 !> \param para_env ...
1095 ! **************************************************************************************************
1096  SUBROUTINE qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, &
1097  nspin, nkp, nkp_local, nmixing, scf_section, para_env)
1098 
1099  TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1100  COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeffs
1101  INTEGER, INTENT(IN) :: ib, nb
1102  REAL(kind=dp), INTENT(IN) :: delta
1103  REAL(kind=dp), INTENT(OUT) :: error_max
1104  LOGICAL, INTENT(OUT) :: diis_step
1105  REAL(kind=dp), INTENT(IN) :: eps_diis
1106  INTEGER, INTENT(IN) :: nspin, nkp, nkp_local
1107  INTEGER, INTENT(IN), OPTIONAL :: nmixing
1108  TYPE(section_vals_type), POINTER :: scf_section
1109  TYPE(mp_para_env_type), POINTER :: para_env
1110 
1111  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_diis_b_step_kp'
1112  REAL(kind=dp), PARAMETER :: eigenvalue_threshold = 1.0e-12_dp
1113 
1114  CHARACTER(LEN=2*default_string_length) :: message
1115  COMPLEX(KIND=dp) :: tmp
1116  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
1117  INTEGER :: handle, ikp, ispin, jb, my_nmixing, nb1, &
1118  output_unit
1119  LOGICAL :: eigenvectors_discarded
1120  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ev
1121  TYPE(cp_cfm_type) :: old_errors
1122  TYPE(cp_cfm_type), POINTER :: new_errors
1123  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1124  TYPE(cp_fm_type) :: ierr, rerr
1125  TYPE(cp_logger_type), POINTER :: logger
1126 
1127  NULLIFY (matrix_struct, new_errors, logger)
1128 
1129  CALL timeset(routinen, handle)
1130 
1131  diis_step = .false.
1132 
1133  my_nmixing = 2
1134  IF (PRESENT(nmixing)) my_nmixing = nmixing
1135 
1136  logger => cp_get_default_logger()
1137 
1138  ! Quick return, if no DIIS is requested
1139  IF (diis_buffer%nbuffer < 1) THEN
1140  CALL timestop(handle)
1141  RETURN
1142  END IF
1143 
1144  ! Check, if a DIIS step is appropriate
1145  diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
1146 
1147  ! Calculate the DIIS buffer, and update it if max_error < eps_diis
1148  CALL cp_cfm_get_info(diis_buffer%error(ib, 1, 1), matrix_struct=matrix_struct)
1149  CALL cp_fm_create(ierr, matrix_struct)
1150  CALL cp_fm_create(rerr, matrix_struct)
1151  CALL cp_cfm_create(old_errors, matrix_struct)
1152  ALLOCATE (b(nb, nb))
1153  b = 0.0_dp
1154  DO jb = 1, nb
1155  DO ikp = 1, nkp_local
1156  DO ispin = 1, nspin
1157  new_errors => diis_buffer%error(ib, ispin, ikp)
1158  CALL cp_cfm_to_fm(diis_buffer%error(jb, ispin, ikp), rerr, ierr)
1159  CALL cp_fm_scale(-1.0_dp, ierr)
1160  CALL cp_fm_to_cfm(rerr, ierr, old_errors)
1161  CALL cp_cfm_trace(old_errors, new_errors, tmp)
1162  b(jb, ib) = b(jb, ib) + 1.0_dp/real(nkp, dp)*tmp
1163  END DO
1164  END DO
1165  b(ib, jb) = conjg(b(jb, ib))
1166  END DO
1167  CALL cp_fm_release(ierr)
1168  CALL cp_fm_release(rerr)
1169  CALL cp_cfm_release(old_errors)
1170  CALL para_env%sum(b)
1171 
1172  error_max = sqrt(real(b(ib, ib))**2 + aimag(b(ib, ib))**2)
1173 
1174  output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
1175  extension=".scfLog")
1176  IF (output_unit > 0) THEN
1177  WRITE (unit=output_unit, fmt="(/,T9,A,I4,/,(T9,A,ES12.3))") &
1178  "DIIS | Current SCF DIIS buffer size: ", nb, &
1179  "DIIS | Maximum SCF DIIS error at last step: ", error_max, &
1180  "DIIS | Current SCF convergence: ", delta, &
1181  "DIIS | Threshold value for a DIIS step: ", eps_diis
1182  IF (error_max < eps_diis) THEN
1183  WRITE (unit=output_unit, fmt="(T9,A)") &
1184  "DIIS | => The SCF DIIS buffer will be updated"
1185  ELSE
1186  WRITE (unit=output_unit, fmt="(T9,A)") &
1187  "DIIS | => No update of the SCF DIIS buffer"
1188  END IF
1189  IF (diis_step .AND. (error_max < eps_diis)) THEN
1190  WRITE (unit=output_unit, fmt="(T9,A,/)") &
1191  "DIIS | => A SCF DIIS step will be performed"
1192  ELSE
1193  WRITE (unit=output_unit, fmt="(T9,A,/)") &
1194  "DIIS | => No SCF DIIS step will be performed"
1195  END IF
1196  END IF
1197 
1198  ! Update the SCF DIIS buffer
1199  IF (error_max < eps_diis) THEN
1200  DO jb = 1, nb
1201  diis_buffer%b_matrix(ib, jb) = b(ib, jb)
1202  diis_buffer%b_matrix(jb, ib) = b(jb, ib)
1203  END DO
1204  ELSE
1205 
1206  diis_step = .false.
1207  END IF
1208  DEALLOCATE (b)
1209 
1210  ! Perform DIIS step
1211  IF (diis_step) THEN
1212 
1213  nb1 = nb + 1
1214 
1215  ALLOCATE (a(nb1, nb1))
1216  ALLOCATE (b(nb1, nb1))
1217  ALLOCATE (ev(nb1))
1218 
1219  ! Set up the linear DIIS equation system
1220  b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
1221 
1222  b(1:nb, nb1) = -1.0_dp
1223  b(nb1, 1:nb) = -1.0_dp
1224  b(nb1, nb1) = 0.0_dp
1225 
1226  ! Solve the linear DIIS equation system
1227  ev(1:nb1) = 0.0_dp !eigenvalues
1228  a(1:nb1, 1:nb1) = 0.0_dp !eigenvectors
1229  CALL complex_diag(b(1:nb1, 1:nb1), a(1:nb1, 1:nb1), ev(1:nb1))
1230  b(1:nb1, 1:nb1) = a(1:nb1, 1:nb1)
1231 
1232  eigenvectors_discarded = .false.
1233 
1234  DO jb = 1, nb1
1235  IF (abs(ev(jb)) < eigenvalue_threshold) THEN
1236  IF (output_unit > 0) THEN
1237  IF (.NOT. eigenvectors_discarded) THEN
1238  WRITE (unit=output_unit, fmt="(T9,A)") &
1239  "DIIS | Checking eigenvalues of the DIIS error matrix"
1240  END IF
1241  WRITE (unit=message, fmt="(T9,A,I6,A,ES10.1,A,ES10.1)") &
1242  "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
1243  "threshold ", eigenvalue_threshold
1244  CALL compress(message)
1245  WRITE (unit=output_unit, fmt="(T9,A)") trim(message)
1246  eigenvectors_discarded = .true.
1247  END IF
1248  a(1:nb1, jb) = 0.0_dp
1249  ELSE
1250  a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
1251  END IF
1252  END DO
1253 
1254  IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
1255  WRITE (unit=output_unit, fmt="(T9,A,/)") &
1256  "DIIS | The corresponding eigenvectors were discarded"
1257  END IF
1258 
1259  coeffs(1:nb) = -matmul(a(1:nb, 1:nb1), conjg(b(nb1, 1:nb1)))
1260  ELSE
1261 
1262  coeffs(:) = 0.0_dp
1263  coeffs(ib) = 1.0_dp
1264  END IF
1265 
1266  CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1267  "PRINT%DIIS_INFO")
1268 
1269  CALL timestop(handle)
1270 
1271  END SUBROUTINE qs_diis_b_step_kp
1272 END MODULE qs_diis
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition: cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_trace
subroutine, public cp_cfm_trace(matrix_a, matrix_b, trace)
Returns the trace of matrix_a^T matrix_b, i.e sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
Definition: cp_cfm_basic_linalg.F:1037
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition: cp_cfm_types.F:817
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition: cp_cfm_types.F:765
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1680
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition: cp_fm_basic_linalg.F:576
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_maxabsval
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(matrix))
Definition: cp_fm_types.F:1064
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
dm_ls_scf_types
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Definition: dm_ls_scf_types.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::complex_diag
subroutine, public complex_diag(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition: mathlib.F:1746
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:372
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
qs_diis_types
buffer for the diis of the scf
Definition: qs_diis_types.F:14
qs_diis
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition: qs_diis.F:21
qs_diis::qs_diis_b_info_kp
subroutine, public qs_diis_b_info_kp(diis_buffer, ib, nb)
Update info about the current buffer step ib and the current number of buffers nb.
Definition: qs_diis.F:987
qs_diis::qs_diis_b_create_sparse
pure subroutine, public qs_diis_b_create_sparse(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer for LS-SCF calculation.
Definition: qs_diis.F:847
qs_diis::qs_diis_b_clear
pure subroutine, public qs_diis_b_clear(diis_buffer)
clears the buffer
Definition: qs_diis.F:517
qs_diis::qs_diis_b_create
subroutine, public qs_diis_b_create(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer.
Definition: qs_diis.F:101
qs_diis::qs_diis_b_step_4lscf
subroutine, public qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, diis_step, eps_diis, nmixing, s_matrix, threshold)
Update the SCF DIIS buffer in linear scaling SCF (LS-SCF), and if appropriate does a diis step.
Definition: qs_diis.F:544
qs_diis::qs_diis_b_calc_err_kp
subroutine, public qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, nkp_local, scf_section)
Calculate and store the error for a given k-point.
Definition: qs_diis.F:1011
qs_diis::qs_diis_b_create_kp
subroutine, public qs_diis_b_create_kp(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer for k-points.
Definition: qs_diis.F:865
qs_diis::qs_diis_b_step_kp
subroutine, public qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, nspin, nkp, nkp_local, nmixing, scf_section, para_env)
Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points.
Definition: qs_diis.F:1098
qs_diis::qs_diis_b_step
subroutine, public qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
Update the SCF DIIS buffer, and if appropriate does a diis step.
Definition: qs_diis.F:228
qs_diis::qs_diis_b_clear_sparse
pure subroutine, public qs_diis_b_clear_sparse(diis_buffer)
clears the DIIS buffer in LS-SCF calculation
Definition: qs_diis.F:831
qs_diis::qs_diis_b_clear_kp
pure subroutine, public qs_diis_b_clear_kp(diis_buffer)
clears the buffer
Definition: qs_diis.F:973
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition: string_utilities.F:3190