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