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