(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
21MODULE 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
90CONTAINS
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
1272END 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_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:55
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:74
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_diis_types::qs_diis_buffer_type_kp
build arrau of pointers to diis buffers in the k-point (complex full matrices) case
Definition qs_diis_types.F:88
qs_diis_types::qs_diis_buffer_type_sparse
build array of pointers to diis buffers for sparse matrix case
Definition qs_diis_types.F:69
qs_diis_types::qs_diis_buffer_type
keeps a buffer with the previous values of s,p,k
Definition qs_diis_types.F:43
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46