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 4532 : 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 4532 : CALL timeset(routineN, handle)
117 :
118 4532 : NULLIFY (diis_buffer%b_matrix)
119 4532 : NULLIFY (diis_buffer%error)
120 4532 : NULLIFY (diis_buffer%param)
121 4532 : diis_buffer%nbuffer = nbuffer
122 4532 : diis_buffer%ncall = 0
123 :
124 4532 : CALL timestop(handle)
125 :
126 4532 : 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 98265 : 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 98265 : CALL timeset(routineN, handle)
160 :
161 98265 : logger => cp_get_default_logger()
162 :
163 98265 : nbuffer = diis_buffer%nbuffer
164 :
165 98265 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
166 34648 : ALLOCATE (diis_buffer%error(nbuffer, nspin))
167 :
168 7620 : DO ispin = 1, nspin
169 24292 : 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 20840 : matrix_struct=matrix_struct)
175 : END DO
176 : END DO
177 : END IF
178 :
179 98265 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
180 34648 : ALLOCATE (diis_buffer%param(nbuffer, nspin))
181 :
182 7620 : DO ispin = 1, nspin
183 24292 : 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 20840 : matrix_struct=matrix_struct)
189 : END DO
190 : END DO
191 : END IF
192 :
193 98265 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
194 13808 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
195 107012 : diis_buffer%b_matrix = 0.0_dp
196 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
197 3452 : extension=".scfLog")
198 3452 : 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 3452 : "PRINT%DIIS_INFO")
204 : END IF
205 :
206 98265 : CALL timestop(handle)
207 :
208 98265 : 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 98265 : 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 98265 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev, occ
258 98265 : REAL(KIND=dp), DIMENSION(:), POINTER :: occa, occb
259 98265 : 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, old_errors, parameters
262 : TYPE(cp_logger_type), POINTER :: logger
263 :
264 : ! -------------------------------------------------------------------------
265 :
266 98265 : CALL timeset(routineN, handle)
267 :
268 98265 : nspin = SIZE(mo_array)
269 98265 : diis_step = .FALSE.
270 :
271 98265 : IF (PRESENT(roks)) THEN
272 1018 : my_roks = .TRUE.
273 1018 : nspin = 1
274 : ELSE
275 : my_roks = .FALSE.
276 : END IF
277 :
278 98265 : my_nmixing = 2
279 98265 : IF (PRESENT(nmixing)) my_nmixing = nmixing
280 :
281 98265 : NULLIFY (c, new_errors, old_errors, parameters, matrix_struct, a, b, occa, occb)
282 98265 : logger => cp_get_default_logger()
283 :
284 : ! Quick return, if no DIIS is requested
285 :
286 98265 : IF (diis_buffer%nbuffer < 1) THEN
287 0 : CALL timestop(handle)
288 : RETURN
289 : END IF
290 :
291 : CALL cp_fm_get_info(kc(1), &
292 98265 : matrix_struct=matrix_struct)
293 : CALL qs_diis_b_check_i_alloc(diis_buffer, &
294 : matrix_struct=matrix_struct, &
295 : nspin=nspin, &
296 98265 : scf_section=scf_section)
297 :
298 98265 : error_max = 0.0_dp
299 :
300 98265 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
301 98265 : diis_buffer%ncall = diis_buffer%ncall + 1
302 98265 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
303 :
304 212288 : DO ispin = 1, nspin
305 :
306 : CALL get_mo_set(mo_set=mo_array(ispin), &
307 : nao=nao, &
308 : nmo=nmo, &
309 : homo=homo, &
310 : mo_coeff=c, &
311 : occupation_numbers=occa, &
312 114023 : maxocc=maxocc)
313 :
314 114023 : new_errors => diis_buffer%error(ib, ispin)
315 114023 : parameters => diis_buffer%param(ib, ispin)
316 :
317 : ! Copy the Kohn-Sham matrix K to the DIIS buffer
318 :
319 114023 : CALL cp_fm_to_fm(kc(ispin), parameters)
320 :
321 114023 : IF (my_roks) THEN
322 :
323 3054 : ALLOCATE (occ(nmo))
324 :
325 : CALL get_mo_set(mo_set=mo_array(2), &
326 1018 : occupation_numbers=occb)
327 :
328 16350 : DO imo = 1, nmo
329 16350 : occ(imo) = SQRT(occa(imo) + occb(imo))
330 : END DO
331 :
332 1018 : CALL cp_fm_to_fm(c, sc)
333 1018 : CALL cp_fm_column_scale(sc, occ(1:homo))
334 :
335 : ! KC <- K*C
336 1018 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, parameters, sc, 0.0_dp, kc(ispin))
337 :
338 1018 : IF (PRESENT(s_matrix)) THEN
339 558 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
340 : ! SC <- S*C
341 558 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, new_errors, c, 0.0_dp, sc)
342 558 : CALL cp_fm_column_scale(sc, occ(1:homo))
343 : END IF
344 :
345 : ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
346 : ! or for an orthogonal basis
347 : ! new_errors <- KC*C^T - C*(KC)^T = K*P - P*K with S = I
348 1018 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
349 1018 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
350 :
351 1018 : DEALLOCATE (occ)
352 :
353 : ELSE
354 :
355 : ! KC <- K*C
356 113005 : CALL cp_fm_symm("L", "U", nao, homo, maxocc, parameters, c, 0.0_dp, kc(ispin))
357 :
358 113005 : IF (PRESENT(s_matrix)) THEN
359 : ! I guess that this copy can be avoided for LSD
360 97029 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
361 : ! sc <- S*C
362 97029 : CALL cp_fm_symm("L", "U", nao, homo, 2.0_dp, new_errors, c, 0.0_dp, sc)
363 : ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
364 97029 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
365 97029 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
366 : ELSE
367 : ! new_errors <- KC*(C)^T - C*(KC)^T = K*P - P*K
368 15976 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, c, kc(ispin), 0.0_dp, new_errors)
369 15976 : CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), c, -1.0_dp, new_errors)
370 : END IF
371 :
372 : END IF
373 :
374 114023 : CALL cp_fm_maxabsval(new_errors, tmp)
375 326311 : error_max = MAX(error_max, tmp)
376 :
377 : END DO
378 :
379 : ! Check, if a DIIS step is appropriate
380 :
381 98265 : diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
382 :
383 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
384 98265 : extension=".scfLog")
385 98265 : IF (output_unit > 0) THEN
386 : WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
387 251 : "DIIS | Current SCF DIIS buffer size: ", nb, &
388 251 : "DIIS | Maximum SCF DIIS error vector element:", error_max, &
389 251 : "DIIS | Current SCF convergence: ", delta, &
390 502 : "DIIS | Threshold value for a DIIS step: ", eps_diis
391 251 : IF (error_max < eps_diis) THEN
392 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
393 102 : "DIIS | => The SCF DIIS buffer will be updated"
394 : ELSE
395 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
396 149 : "DIIS | => No update of the SCF DIIS buffer"
397 : END IF
398 251 : IF (diis_step .AND. (error_max < eps_diis)) THEN
399 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
400 62 : "DIIS | => A SCF DIIS step will be performed"
401 : ELSE
402 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
403 189 : "DIIS | => No SCF DIIS step will be performed"
404 : END IF
405 : END IF
406 :
407 : ! Update the SCF DIIS buffer
408 :
409 98265 : IF (error_max < eps_diis) THEN
410 :
411 84577 : b => diis_buffer%b_matrix
412 :
413 358895 : DO jb = 1, nb
414 274318 : b(jb, ib) = 0.0_dp
415 596230 : DO ispin = 1, nspin
416 321912 : old_errors => diis_buffer%error(jb, ispin)
417 321912 : new_errors => diis_buffer%error(ib, ispin)
418 321912 : CALL cp_fm_trace(old_errors, new_errors, tmp)
419 596230 : b(jb, ib) = b(jb, ib) + tmp
420 : END DO
421 358895 : b(ib, jb) = b(jb, ib)
422 : END DO
423 :
424 : ELSE
425 :
426 13688 : diis_step = .FALSE.
427 :
428 : END IF
429 :
430 : ! Perform DIIS step
431 :
432 98265 : IF (diis_step) THEN
433 :
434 62895 : nb1 = nb + 1
435 :
436 251580 : ALLOCATE (a(nb1, nb1))
437 188685 : ALLOCATE (b(nb1, nb1))
438 188685 : ALLOCATE (ev(nb1))
439 :
440 : ! Set up the linear DIIS equation system
441 :
442 2215279 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
443 :
444 289711 : b(1:nb, nb1) = -1.0_dp
445 289711 : b(nb1, 1:nb) = -1.0_dp
446 62895 : b(nb1, nb1) = 0.0_dp
447 :
448 : ! Solve the linear DIIS equation system
449 :
450 352606 : ev(1:nb1) = 0.0_dp
451 62895 : CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
452 :
453 3374123 : a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
454 :
455 62895 : eigenvectors_discarded = .FALSE.
456 :
457 352606 : DO jb = 1, nb1
458 352606 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
459 44256 : IF (output_unit > 0) THEN
460 5 : IF (.NOT. eigenvectors_discarded) THEN
461 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
462 5 : "DIIS | Checking eigenvalues of the DIIS error matrix"
463 : END IF
464 : WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
465 5 : "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
466 10 : "threshold ", eigenvalue_threshold
467 5 : CALL compress(message)
468 5 : WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
469 5 : eigenvectors_discarded = .TRUE.
470 : END IF
471 259036 : a(1:nb1, jb) = 0.0_dp
472 : ELSE
473 1396578 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
474 : END IF
475 : END DO
476 :
477 62895 : IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
478 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
479 5 : "DIIS | The corresponding eigenvectors were discarded"
480 : END IF
481 :
482 1655614 : ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
483 :
484 : ! Update Kohn-Sham matrix
485 :
486 136308 : DO ispin = 1, nspin
487 73413 : CALL cp_fm_set_all(kc(ispin), 0.0_dp)
488 402546 : DO jb = 1, nb
489 266238 : parameters => diis_buffer%param(jb, ispin)
490 339651 : CALL cp_fm_scale_and_add(1.0_dp, kc(ispin), -ev(jb), parameters)
491 : END DO
492 : END DO
493 :
494 62895 : DEALLOCATE (a)
495 62895 : DEALLOCATE (b)
496 62895 : DEALLOCATE (ev)
497 :
498 : ELSE
499 :
500 75980 : DO ispin = 1, nspin
501 40610 : parameters => diis_buffer%param(ib, ispin)
502 75980 : CALL cp_fm_to_fm(parameters, kc(ispin))
503 : END DO
504 :
505 : END IF
506 :
507 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
508 98265 : "PRINT%DIIS_INFO")
509 :
510 98265 : CALL timestop(handle)
511 :
512 196530 : END SUBROUTINE qs_diis_b_step
513 :
514 : ! **************************************************************************************************
515 : !> \brief clears the buffer
516 : !> \param diis_buffer the buffer to clear
517 : !> \par History
518 : !> 02.2003 created [fawzi]
519 : !> \author fawzi
520 : ! **************************************************************************************************
521 16410 : PURE SUBROUTINE qs_diis_b_clear(diis_buffer)
522 :
523 : TYPE(qs_diis_buffer_type), INTENT(INOUT) :: diis_buffer
524 :
525 16410 : diis_buffer%ncall = 0
526 :
527 16410 : END SUBROUTINE qs_diis_b_clear
528 :
529 : ! **************************************************************************************************
530 : !> \brief Update the SCF DIIS buffer in linear scaling SCF (LS-SCF),
531 : !> and if appropriate does a diis step.
532 : !> \param diis_buffer ...
533 : !> \param qs_env ...
534 : !> \param ls_scf_env ...
535 : !> \param unit_nr ...
536 : !> \param iscf ...
537 : !> \param diis_step ...
538 : !> \param eps_diis ...
539 : !> \param nmixing ...
540 : !> \param s_matrix ...
541 : !> \param threshold ...
542 : !> \par History
543 : !> - Adapted for LS-SCF (10-11-14) from qs_diis_b_step
544 : !> \author Fredy W. Aquino
545 : ! **************************************************************************************************
546 :
547 18 : SUBROUTINE qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, &
548 : diis_step, eps_diis, nmixing, s_matrix, threshold)
549 : ! Note.- Input: ls_scf_env%matrix_p(ispin) , Density Matrix
550 : ! matrix_ks (from qs_env) , Kohn-Sham Matrix (IN/OUT)
551 :
552 : TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer
553 : TYPE(qs_environment_type), POINTER :: qs_env
554 : TYPE(ls_scf_env_type) :: ls_scf_env
555 : INTEGER, INTENT(IN) :: unit_nr, iscf
556 : LOGICAL, INTENT(OUT) :: diis_step
557 : REAL(KIND=dp), INTENT(IN) :: eps_diis
558 : INTEGER, INTENT(IN), OPTIONAL :: nmixing
559 : TYPE(dbcsr_type), OPTIONAL :: s_matrix
560 : REAL(KIND=dp), INTENT(IN) :: threshold
561 :
562 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_4lscf'
563 : REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
564 :
565 : INTEGER :: handle, ib, ispin, jb, my_nmixing, nb, &
566 : nb1, nspin
567 : REAL(KIND=dp) :: error_max, tmp
568 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
569 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
570 : TYPE(cp_logger_type), POINTER :: logger
571 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
572 : TYPE(dbcsr_type) :: matrix_KSerr_t, matrix_tmp
573 : TYPE(dbcsr_type), POINTER :: new_errors, old_errors, parameters
574 : TYPE(mp_para_env_type), POINTER :: para_env
575 :
576 18 : CALL timeset(routineN, handle)
577 18 : IF (ls_scf_env%do_pao) THEN
578 0 : CPABORT("LS_SCF%LS_DIIS not compatible with PAO")
579 : END IF
580 18 : nspin = ls_scf_env%nspins
581 18 : diis_step = .FALSE.
582 18 : my_nmixing = 2
583 18 : IF (PRESENT(nmixing)) my_nmixing = nmixing
584 18 : NULLIFY (new_errors, old_errors, parameters, a, b)
585 18 : logger => cp_get_default_logger()
586 : ! Quick return, if no DIIS is requested
587 18 : IF (diis_buffer%nbuffer < 1) THEN
588 0 : CALL timestop(handle)
589 : RETURN
590 : END IF
591 :
592 : ! Getting current Kohn-Sham matrix from qs_env
593 : CALL get_qs_env(qs_env, &
594 : para_env=para_env, &
595 18 : matrix_ks=matrix_ks)
596 : CALL qs_diis_b_check_i_alloc_sparse( &
597 : diis_buffer, &
598 : ls_scf_env, &
599 18 : nspin)
600 18 : error_max = 0.0_dp
601 :
602 18 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
603 18 : diis_buffer%ncall = diis_buffer%ncall + 1
604 18 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
605 : ! Create scratch arrays
606 : CALL dbcsr_create(matrix_tmp, &
607 : template=ls_scf_env%matrix_ks(1), &
608 18 : matrix_type='N')
609 18 : CALL dbcsr_set(matrix_tmp, 0.0_dp) ! reset matrix
610 : CALL dbcsr_create(matrix_KSerr_t, &
611 : template=ls_scf_env%matrix_ks(1), &
612 18 : matrix_type='N')
613 18 : CALL dbcsr_set(matrix_KSerr_t, 0.0_dp) ! reset matrix
614 :
615 46 : DO ispin = 1, nspin ! ------ Loop-ispin----START
616 :
617 28 : new_errors => diis_buffer%error(ib, ispin)%matrix
618 28 : parameters => diis_buffer%param(ib, ispin)%matrix
619 : ! Copy the Kohn-Sham matrix K to the DIIS buffer
620 : CALL dbcsr_copy(parameters, & ! out
621 28 : matrix_ks(ispin)%matrix) ! in
622 :
623 28 : IF (PRESENT(s_matrix)) THEN ! if-s_matrix ---------- START
624 : ! Calculate Kohn-Sham error (non-orthogonal)= K*P*S-(K*P*S)^T
625 : ! matrix_tmp = P*S
626 : CALL dbcsr_multiply("N", "N", &
627 : 1.0_dp, ls_scf_env%matrix_p(ispin), &
628 : s_matrix, &
629 : 0.0_dp, matrix_tmp, &
630 28 : filter_eps=threshold)
631 : ! new_errors= K*P*S
632 : CALL dbcsr_multiply("N", "N", &
633 : 1.0_dp, matrix_ks(ispin)%matrix, &
634 : matrix_tmp, &
635 : 0.0_dp, new_errors, &
636 28 : filter_eps=threshold)
637 : ! matrix_KSerr_t= transpose(K*P*S)
638 : CALL dbcsr_transposed(matrix_KSerr_t, &
639 28 : new_errors)
640 : ! new_errors=K*P*S-transpose(K*P*S)
641 : CALL dbcsr_add(new_errors, &
642 : matrix_KSerr_t, &
643 28 : 1.0_dp, -1.0_dp)
644 : ELSE ! if-s_matrix ---------- MID
645 : ! Calculate Kohn-Sham error (orthogonal)= K*P - P*K
646 : ! new_errors=K*P
647 : CALL dbcsr_multiply("N", "N", &
648 : 1.0_dp, matrix_ks(ispin)%matrix, &
649 : ls_scf_env%matrix_p(ispin), &
650 : 0.0_dp, new_errors, &
651 0 : filter_eps=threshold)
652 : ! matrix_KSerr_t= transpose(K*P)
653 : CALL dbcsr_transposed(matrix_KSerr_t, &
654 0 : new_errors)
655 : ! new_errors=K*P-transpose(K*P)
656 : CALL dbcsr_add(new_errors, &
657 : matrix_KSerr_t, &
658 0 : 1.0_dp, -1.0_dp)
659 : END IF ! if-s_matrix ---------- END
660 :
661 28 : tmp = dbcsr_maxabs(new_errors)
662 46 : error_max = MAX(error_max, tmp)
663 :
664 : END DO ! ------ Loop-ispin----END
665 :
666 : ! Check, if a DIIS step is appropriate
667 :
668 18 : diis_step = (diis_buffer%ncall >= my_nmixing)
669 :
670 18 : IF (unit_nr > 0) THEN
671 : WRITE (unit_nr, '(A29,I3,A3,4(I3,A1))') &
672 9 : "DIIS: (ncall,nbuffer,ib,nb)=(", iscf, ")=(", &
673 18 : diis_buffer%ncall, ",", diis_buffer%nbuffer, ",", ib, ",", nb, ")"
674 : WRITE (unit_nr, '(A57,I3,A3,L1,A1,F10.8,A1,F4.2,A1,L1,A1)') &
675 9 : "DIIS: (diis_step,error_max,eps_diis,error_max<eps_diis)=(", &
676 9 : iscf, ")=(", diis_step, ",", error_max, ",", eps_diis, ",", &
677 18 : (error_max < eps_diis), ")"
678 : WRITE (unit_nr, '(A75)') &
679 9 : "DIIS: diis_step=T : Perform DIIS error_max<eps_diis=T : Update DIIS buffer"
680 : END IF
681 :
682 : ! Update the SCF DIIS buffer
683 18 : IF (error_max < eps_diis) THEN
684 18 : b => diis_buffer%b_matrix
685 66 : DO jb = 1, nb
686 48 : b(jb, ib) = 0.0_dp
687 124 : DO ispin = 1, nspin
688 76 : old_errors => diis_buffer%error(jb, ispin)%matrix
689 76 : new_errors => diis_buffer%error(ib, ispin)%matrix
690 : CALL dbcsr_dot(old_errors, &
691 : new_errors, &
692 76 : tmp) ! out : < f_i | f_j >
693 124 : b(jb, ib) = b(jb, ib) + tmp
694 : END DO ! end-loop-ispin
695 66 : b(ib, jb) = b(jb, ib)
696 : END DO ! end-loop-jb
697 : ELSE
698 0 : diis_step = .FALSE.
699 : END IF
700 :
701 : ! Perform DIIS step
702 18 : IF (diis_step) THEN
703 14 : nb1 = nb + 1
704 56 : ALLOCATE (a(nb1, nb1))
705 42 : ALLOCATE (b(nb1, nb1))
706 42 : ALLOCATE (ev(nb1))
707 : ! Set up the linear DIIS equation system
708 398 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
709 58 : b(1:nb, nb1) = -1.0_dp
710 58 : b(nb1, 1:nb) = -1.0_dp
711 14 : b(nb1, nb1) = 0.0_dp
712 : ! Solve the linear DIIS equation system
713 14 : CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
714 630 : a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
715 72 : DO jb = 1, nb1
716 72 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
717 0 : a(1:nb1, jb) = 0.0_dp
718 : ELSE
719 308 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
720 : END IF
721 : END DO ! end-loop-jb
722 :
723 308 : ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
724 :
725 : ! Update Kohn-Sham matrix
726 14 : IF (iscf >= ls_scf_env%iter_ini_diis) THEN ! if-iscf-to-updateKS------ START
727 :
728 14 : IF (unit_nr > 0) THEN
729 7 : WRITE (unit_nr, '(A40,I3)') 'DIIS: Updating Kohn-Sham matrix at iscf=', iscf
730 : END IF
731 :
732 36 : DO ispin = 1, nspin
733 : CALL dbcsr_set(matrix_ks(ispin)%matrix, & ! reset matrix
734 22 : 0.0_dp)
735 106 : DO jb = 1, nb
736 70 : parameters => diis_buffer%param(jb, ispin)%matrix
737 : CALL dbcsr_add(matrix_ks(ispin)%matrix, parameters, &
738 92 : 1.0_dp, -ev(jb))
739 : END DO ! end-loop-jb
740 : END DO ! end-loop-ispin
741 : END IF ! if-iscf-to-updateKS------ END
742 :
743 14 : DEALLOCATE (a)
744 14 : DEALLOCATE (b)
745 14 : DEALLOCATE (ev)
746 :
747 : ELSE
748 10 : DO ispin = 1, nspin
749 6 : parameters => diis_buffer%param(ib, ispin)%matrix
750 : CALL dbcsr_copy(parameters, & ! out
751 10 : matrix_ks(ispin)%matrix) ! in
752 : END DO ! end-loop-ispin
753 : END IF
754 18 : CALL dbcsr_release(matrix_tmp)
755 18 : CALL dbcsr_release(matrix_KSerr_t)
756 18 : CALL timestop(handle)
757 :
758 18 : END SUBROUTINE qs_diis_b_step_4lscf
759 :
760 : ! **************************************************************************************************
761 : !> \brief Allocate and initialize a DIIS buffer with a buffer size of nbuffer.
762 : !> \param diis_buffer the buffer to initialize
763 : !> \param ls_scf_env ...
764 : !> \param nspin ...
765 : !> \par History
766 : !> - Adapted from qs_diis_b_check_i_alloc for sparse matrices and
767 : !> used in LS-SCF module (ls_scf_main) (10-11-14)
768 : !> \author Fredy W. Aquino
769 : !> \note
770 : !> check to allocate matrices only when needed
771 : ! **************************************************************************************************
772 :
773 18 : SUBROUTINE qs_diis_b_check_i_alloc_sparse(diis_buffer, ls_scf_env, &
774 : nspin)
775 :
776 : TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
777 : TYPE(ls_scf_env_type) :: ls_scf_env
778 : INTEGER, INTENT(IN) :: nspin
779 :
780 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_sparse'
781 :
782 : INTEGER :: handle, ibuffer, ispin, nbuffer
783 : TYPE(cp_logger_type), POINTER :: logger
784 :
785 : ! -------------------------------------------------------------------------
786 :
787 18 : CALL timeset(routineN, handle)
788 :
789 18 : logger => cp_get_default_logger()
790 :
791 18 : nbuffer = diis_buffer%nbuffer
792 :
793 18 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
794 46 : ALLOCATE (diis_buffer%error(nbuffer, nspin))
795 :
796 10 : DO ispin = 1, nspin
797 34 : DO ibuffer = 1, nbuffer
798 24 : ALLOCATE (diis_buffer%error(ibuffer, ispin)%matrix)
799 :
800 : CALL dbcsr_create(diis_buffer%error(ibuffer, ispin)%matrix, &
801 : template=ls_scf_env%matrix_ks(1), &
802 30 : matrix_type='N')
803 : END DO
804 : END DO
805 : END IF
806 :
807 18 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
808 46 : ALLOCATE (diis_buffer%param(nbuffer, nspin))
809 :
810 10 : DO ispin = 1, nspin
811 34 : DO ibuffer = 1, nbuffer
812 24 : ALLOCATE (diis_buffer%param(ibuffer, ispin)%matrix)
813 : CALL dbcsr_create(diis_buffer%param(ibuffer, ispin)%matrix, &
814 : template=ls_scf_env%matrix_ks(1), &
815 30 : matrix_type='N')
816 : END DO
817 : END DO
818 : END IF
819 :
820 18 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
821 16 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
822 :
823 124 : diis_buffer%b_matrix = 0.0_dp
824 : END IF
825 :
826 18 : CALL timestop(handle)
827 :
828 18 : END SUBROUTINE qs_diis_b_check_i_alloc_sparse
829 :
830 : ! **************************************************************************************************
831 : !> \brief clears the DIIS buffer in LS-SCF calculation
832 : !> \param diis_buffer the buffer to clear
833 : !> \par History
834 : !> 10-11-14 created [FA] modified from qs_diis_b_clear
835 : !> \author Fredy W. Aquino
836 : ! **************************************************************************************************
837 :
838 4 : PURE SUBROUTINE qs_diis_b_clear_sparse(diis_buffer)
839 :
840 : TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT) :: diis_buffer
841 :
842 4 : diis_buffer%ncall = 0
843 :
844 4 : END SUBROUTINE qs_diis_b_clear_sparse
845 :
846 : ! **************************************************************************************************
847 : !> \brief Allocates an SCF DIIS buffer for LS-SCF calculation
848 : !> \param diis_buffer the buffer to create
849 : !> \param nbuffer ...
850 : !> \par History
851 : !> 10-11-14 created [FA] modified from qs_diis_b_create
852 : !> \author Fredy W. Aquino
853 : ! **************************************************************************************************
854 4 : PURE SUBROUTINE qs_diis_b_create_sparse(diis_buffer, nbuffer)
855 :
856 : TYPE(qs_diis_buffer_type_sparse), INTENT(OUT) :: diis_buffer
857 : INTEGER, INTENT(in) :: nbuffer
858 :
859 : NULLIFY (diis_buffer%b_matrix)
860 : NULLIFY (diis_buffer%error)
861 : NULLIFY (diis_buffer%param)
862 4 : diis_buffer%nbuffer = nbuffer
863 4 : diis_buffer%ncall = 0
864 :
865 4 : END SUBROUTINE qs_diis_b_create_sparse
866 :
867 : ! **************************************************************************************************
868 : !> \brief Allocates an SCF DIIS buffer for k-points
869 : !> \param diis_buffer the buffer to create
870 : !> \param nbuffer ...
871 : ! **************************************************************************************************
872 2560 : SUBROUTINE qs_diis_b_create_kp(diis_buffer, nbuffer)
873 :
874 : TYPE(qs_diis_buffer_type_kp), INTENT(OUT) :: diis_buffer
875 : INTEGER, INTENT(in) :: nbuffer
876 :
877 : NULLIFY (diis_buffer%b_matrix)
878 : NULLIFY (diis_buffer%error)
879 : NULLIFY (diis_buffer%param)
880 : NULLIFY (diis_buffer%smat)
881 2560 : diis_buffer%nbuffer = nbuffer
882 2560 : diis_buffer%ncall = 0
883 :
884 2560 : END SUBROUTINE qs_diis_b_create_kp
885 :
886 : ! **************************************************************************************************
887 : !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
888 : !> variables and with a buffer size of nbuffer, in the k-point case
889 : !> \param diis_buffer the buffer to initialize
890 : !> \param matrix_struct the structure for the matrix of the buffer note: this is in the kp subgroup
891 : !> \param nspin ...
892 : !> \param nkp ...
893 : !> \param scf_section ...
894 : ! **************************************************************************************************
895 27498 : SUBROUTINE qs_diis_b_check_i_alloc_kp(diis_buffer, matrix_struct, nspin, nkp, scf_section)
896 :
897 : TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
898 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
899 : INTEGER, INTENT(IN) :: nspin, nkp
900 : TYPE(section_vals_type), POINTER :: scf_section
901 :
902 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_kp'
903 :
904 : INTEGER :: handle, ibuffer, ikp, ispin, nbuffer, &
905 : output_unit
906 : TYPE(cp_logger_type), POINTER :: logger
907 :
908 : ! -------------------------------------------------------------------------
909 :
910 27498 : CALL timeset(routineN, handle)
911 :
912 27498 : logger => cp_get_default_logger()
913 :
914 27498 : nbuffer = diis_buffer%nbuffer
915 :
916 27498 : IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
917 35863 : ALLOCATE (diis_buffer%error(nbuffer, nspin, nkp))
918 :
919 5949 : DO ikp = 1, nkp
920 10271 : DO ispin = 1, nspin
921 25483 : DO ibuffer = 1, nbuffer
922 : CALL cp_cfm_create(diis_buffer%error(ibuffer, ispin, ikp), &
923 : name="qs_diis_b%error("// &
924 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
925 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
926 21610 : matrix_struct=matrix_struct)
927 : END DO
928 : END DO
929 : END DO
930 : END IF
931 :
932 27498 : IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
933 35863 : ALLOCATE (diis_buffer%param(nbuffer, nspin, nkp))
934 :
935 5949 : DO ikp = 1, nkp
936 10271 : DO ispin = 1, nspin
937 25483 : DO ibuffer = 1, nbuffer
938 : CALL cp_cfm_create(diis_buffer%param(ibuffer, ispin, ikp), &
939 : name="qs_diis_b%param("// &
940 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
941 : TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
942 21610 : matrix_struct=matrix_struct)
943 : END DO
944 : END DO
945 : END DO
946 : END IF
947 :
948 27498 : IF (.NOT. ASSOCIATED(diis_buffer%smat)) THEN
949 10101 : ALLOCATE (diis_buffer%smat(nkp))
950 5949 : DO ikp = 1, nkp
951 : CALL cp_cfm_create(diis_buffer%smat(ikp), &
952 : name="kp_cfm_smat("//TRIM(cp_to_string(ikp))//")", &
953 5949 : matrix_struct=matrix_struct)
954 : END DO
955 : END IF
956 :
957 27498 : IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
958 8304 : ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
959 64356 : diis_buffer%b_matrix = 0.0_dp
960 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
961 2076 : extension=".scfLog")
962 2076 : IF (output_unit > 0) THEN
963 : WRITE (UNIT=output_unit, FMT="(/,T9,A)") &
964 0 : "DIIS | The SCF DIIS buffer was allocated and initialized"
965 : END IF
966 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
967 2076 : "PRINT%DIIS_INFO")
968 : END IF
969 :
970 27498 : CALL timestop(handle)
971 :
972 27498 : END SUBROUTINE qs_diis_b_check_i_alloc_kp
973 :
974 : ! **************************************************************************************************
975 : !> \brief clears the buffer
976 : !> \param diis_buffer the buffer to clear
977 : ! **************************************************************************************************
978 3210 : PURE SUBROUTINE qs_diis_b_clear_kp(diis_buffer)
979 :
980 : TYPE(qs_diis_buffer_type_kp), INTENT(INOUT) :: diis_buffer
981 :
982 3210 : diis_buffer%ncall = 0
983 :
984 3210 : END SUBROUTINE qs_diis_b_clear_kp
985 :
986 : ! **************************************************************************************************
987 : !> \brief Update info about the current buffer step ib and the current number of buffers nb
988 : !> \param diis_buffer ...
989 : !> \param ib ...
990 : !> \param nb ...
991 : ! **************************************************************************************************
992 27498 : SUBROUTINE qs_diis_b_info_kp(diis_buffer, ib, nb)
993 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
994 : INTEGER, INTENT(OUT) :: ib, nb
995 :
996 27498 : ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
997 27498 : diis_buffer%ncall = diis_buffer%ncall + 1
998 27498 : nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
999 :
1000 27498 : END SUBROUTINE qs_diis_b_info_kp
1001 :
1002 : ! **************************************************************************************************
1003 : !> \brief Calculate and store the error for a given k-point
1004 : !> \param diis_buffer ...
1005 : !> \param ib ...
1006 : !> \param mos ...
1007 : !> \param kc ...
1008 : !> \param sc ...
1009 : !> \param ispin ...
1010 : !> \param ikp ...
1011 : !> \param cmos caller-owned scratch; the controller preallocates it and the DIIS buffer
1012 : !> \param real_wfn ...
1013 : !> \param local use local GEMM on a one-rank k-point group
1014 : !> \param gemm_ctx private local GEMM context, unused for distributed groups
1015 : !> \note We assume that we always have an overlap matrix.
1016 : !> TODO: do we need to pass the kp weight for the back Fourier transform?
1017 : ! **************************************************************************************************
1018 97664 : SUBROUTINE qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, cmos, real_wfn, local, gemm_ctx)
1019 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1020 : INTEGER, INTENT(IN) :: ib
1021 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
1022 : TYPE(cp_cfm_type), INTENT(INOUT) :: kc, sc
1023 : INTEGER, INTENT(IN) :: ispin, ikp
1024 : TYPE(cp_cfm_type), INTENT(INOUT) :: cmos
1025 : LOGICAL, INTENT(IN), OPTIONAL :: real_wfn, local
1026 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: gemm_ctx
1027 :
1028 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_calc_err_kp'
1029 :
1030 : INTEGER :: handle, homo, nao
1031 : LOGICAL :: my_real_wfn
1032 : REAL(dp) :: maxocc
1033 : TYPE(cp_cfm_type), POINTER :: new_errors, parameters, smat
1034 : TYPE(cp_fm_type), POINTER :: imos, rmos
1035 :
1036 48832 : NULLIFY (imos, rmos, parameters, new_errors, smat)
1037 :
1038 48832 : CALL timeset(routineN, handle)
1039 48832 : my_real_wfn = .FALSE.
1040 48832 : IF (PRESENT(real_wfn)) my_real_wfn = real_wfn
1041 :
1042 : !Calculate the error for this given k-point, store the KS matrix as well as the ovlp matrix
1043 : !All of this happens within the kp subgroups
1044 :
1045 : ! Quick return, if no DIIS is requested
1046 48832 : IF (diis_buffer%nbuffer < 1) THEN
1047 0 : CALL timestop(handle)
1048 0 : RETURN
1049 : END IF
1050 48832 : CALL get_mo_set(mos(1, ispin), nao=nao, homo=homo, mo_coeff=rmos, maxocc=maxocc)
1051 :
1052 48832 : IF (my_real_wfn) THEN
1053 : CALL cp_cfm_scale_and_add_fm(CMPLX(0.0_dp, KIND=dp), cmos, &
1054 0 : CMPLX(1.0_dp, KIND=dp), rmos)
1055 : ELSE
1056 48832 : CALL get_mo_set(mos(2, ispin), mo_coeff=imos)
1057 48832 : CALL cp_fm_to_cfm(rmos, imos, cmos)
1058 : END IF
1059 :
1060 48832 : new_errors => diis_buffer%error(ib, ispin, ikp)
1061 48832 : parameters => diis_buffer%param(ib, ispin, ikp)
1062 48832 : smat => diis_buffer%smat(ikp)
1063 :
1064 : !copy the KS and overlap matrices to the DIIS buffer
1065 48832 : CALL cp_cfm_to_cfm(kc, parameters)
1066 48832 : CALL cp_cfm_to_cfm(sc, smat)
1067 :
1068 : ! KC <- K*C; SC <- S*C. Retain the existing factor of two in SC.
1069 : ! new_errors <- KC*(SC)^H - SC*(KC)^H.
1070 48832 : IF (local) THEN
1071 100734 : CPASSERT(PRODUCT(parameters%matrix_struct%context%num_pe) == 1)
1072 100734 : CPASSERT(PRODUCT(cmos%matrix_struct%context%num_pe) == 1)
1073 100734 : CPASSERT(PRODUCT(kc%matrix_struct%context%num_pe) == 1)
1074 100734 : CPASSERT(PRODUCT(smat%matrix_struct%context%num_pe) == 1)
1075 100734 : CPASSERT(PRODUCT(sc%matrix_struct%context%num_pe) == 1)
1076 100734 : CPASSERT(PRODUCT(new_errors%matrix_struct%context%num_pe) == 1)
1077 : CALL gemm_ctx%gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), &
1078 : parameters%local_data, SIZE(parameters%local_data, 1), &
1079 : cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
1080 33578 : kc%local_data, SIZE(kc%local_data, 1))
1081 : CALL gemm_ctx%gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), &
1082 : smat%local_data, SIZE(smat%local_data, 1), &
1083 : cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
1084 33578 : sc%local_data, SIZE(sc%local_data, 1))
1085 : CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
1086 : sc%local_data, SIZE(sc%local_data, 1), &
1087 : kc%local_data, SIZE(kc%local_data, 1), (0.0_dp, 0.0_dp), &
1088 33578 : new_errors%local_data, SIZE(new_errors%local_data, 1))
1089 : CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
1090 : kc%local_data, SIZE(kc%local_data, 1), &
1091 : sc%local_data, SIZE(sc%local_data, 1), (-1.0_dp, 0.0_dp), &
1092 33578 : new_errors%local_data, SIZE(new_errors%local_data, 1))
1093 : ELSE
1094 15254 : CALL parallel_gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), parameters, cmos, (0.0_dp, 0.0_dp), kc)
1095 15254 : CALL parallel_gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), smat, cmos, (0.0_dp, 0.0_dp), sc)
1096 15254 : CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), sc, kc, (0.0_dp, 0.0_dp), new_errors)
1097 15254 : CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), kc, sc, (-1.0_dp, 0.0_dp), new_errors)
1098 : END IF
1099 :
1100 48832 : CALL timestop(handle)
1101 :
1102 : END SUBROUTINE qs_diis_b_calc_err_kp
1103 :
1104 : ! **************************************************************************************************
1105 : !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points
1106 : !> \param diis_buffer ...
1107 : !> \param coeffs ...
1108 : !> \param ib ...
1109 : !> \param nb ...
1110 : !> \param delta ...
1111 : !> \param error_max ...
1112 : !> \param diis_step ...
1113 : !> \param eps_diis ...
1114 : !> \param nspin ...
1115 : !> \param nkp ...
1116 : !> \param nkp_local ...
1117 : !> \param nmixing ...
1118 : !> \param scf_section ...
1119 : !> \param para_env_inter_kp communicator connecting the k-point groups
1120 : ! **************************************************************************************************
1121 27498 : SUBROUTINE qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, &
1122 : nspin, nkp, nkp_local, nmixing, scf_section, para_env_inter_kp)
1123 :
1124 : TYPE(qs_diis_buffer_type_kp), POINTER :: diis_buffer
1125 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeffs
1126 : INTEGER, INTENT(IN) :: ib, nb
1127 : REAL(KIND=dp), INTENT(IN) :: delta
1128 : REAL(KIND=dp), INTENT(OUT) :: error_max
1129 : LOGICAL, INTENT(OUT) :: diis_step
1130 : REAL(KIND=dp), INTENT(IN) :: eps_diis
1131 : INTEGER, INTENT(IN) :: nspin, nkp, nkp_local
1132 : INTEGER, INTENT(IN), OPTIONAL :: nmixing
1133 : TYPE(section_vals_type), POINTER :: scf_section
1134 : TYPE(mp_para_env_type), POINTER :: para_env_inter_kp
1135 :
1136 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_kp'
1137 : REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
1138 :
1139 : CHARACTER(LEN=2*default_string_length) :: message
1140 : COMPLEX(KIND=dp) :: tmp
1141 27498 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
1142 : INTEGER :: handle, ikp, ispin, jb, my_nmixing, nb1, &
1143 : output_unit
1144 : LOGICAL :: eigenvectors_discarded
1145 27498 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
1146 : TYPE(cp_cfm_type) :: old_errors
1147 : TYPE(cp_cfm_type), POINTER :: new_errors
1148 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1149 : TYPE(cp_fm_type) :: ierr, rerr
1150 : TYPE(cp_logger_type), POINTER :: logger
1151 :
1152 27498 : NULLIFY (matrix_struct, new_errors, logger)
1153 :
1154 27498 : CALL timeset(routineN, handle)
1155 :
1156 27498 : diis_step = .FALSE.
1157 :
1158 27498 : my_nmixing = 2
1159 27498 : IF (PRESENT(nmixing)) my_nmixing = nmixing
1160 :
1161 27498 : logger => cp_get_default_logger()
1162 :
1163 : ! Quick return, if no DIIS is requested
1164 27498 : IF (diis_buffer%nbuffer < 1) THEN
1165 0 : CALL timestop(handle)
1166 0 : RETURN
1167 : END IF
1168 :
1169 : ! Check, if a DIIS step is appropriate
1170 27498 : diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
1171 :
1172 : ! Calculate the DIIS buffer, and update it if max_error < eps_diis
1173 27498 : CALL cp_cfm_get_info(diis_buffer%error(ib, 1, 1), matrix_struct=matrix_struct)
1174 27498 : CALL cp_fm_create(ierr, matrix_struct)
1175 27498 : CALL cp_fm_create(rerr, matrix_struct)
1176 27498 : CALL cp_cfm_create(old_errors, matrix_struct)
1177 109992 : ALLOCATE (b(nb, nb))
1178 27498 : b = 0.0_dp
1179 122794 : DO jb = 1, nb
1180 243442 : DO ikp = 1, nkp_local
1181 407349 : DO ispin = 1, nspin
1182 163907 : new_errors => diis_buffer%error(ib, ispin, ikp)
1183 163907 : CALL cp_cfm_to_fm(diis_buffer%error(jb, ispin, ikp), rerr, ierr)
1184 163907 : CALL cp_fm_scale(-1.0_dp, ierr)
1185 163907 : CALL cp_fm_to_cfm(rerr, ierr, old_errors)
1186 163907 : CALL cp_cfm_trace(old_errors, new_errors, tmp)
1187 312053 : b(jb, ib) = b(jb, ib) + 1.0_dp/REAL(nkp, dp)*tmp
1188 : END DO
1189 : END DO
1190 122794 : b(ib, jb) = CONJG(b(jb, ib))
1191 : END DO
1192 27498 : CALL cp_fm_release(ierr)
1193 27498 : CALL cp_fm_release(rerr)
1194 27498 : CALL cp_cfm_release(old_errors)
1195 : ! cp_cfm_trace has already reduced each trace over the ranks of its
1196 : ! k-point group. Reduce the group contributions only once, over the
1197 : ! communicator that connects the k-point groups.
1198 27498 : CALL para_env_inter_kp%sum(b)
1199 :
1200 27498 : error_max = SQRT(REAL(b(ib, ib))**2 + AIMAG(b(ib, ib))**2)
1201 :
1202 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
1203 27498 : extension=".scfLog")
1204 27498 : IF (output_unit > 0) THEN
1205 : WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
1206 0 : "DIIS | Current SCF DIIS buffer size: ", nb, &
1207 0 : "DIIS | Maximum SCF DIIS error at last step: ", error_max, &
1208 0 : "DIIS | Current SCF convergence: ", delta, &
1209 0 : "DIIS | Threshold value for a DIIS step: ", eps_diis
1210 0 : IF (error_max < eps_diis) THEN
1211 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1212 0 : "DIIS | => The SCF DIIS buffer will be updated"
1213 : ELSE
1214 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1215 0 : "DIIS | => No update of the SCF DIIS buffer"
1216 : END IF
1217 0 : IF (diis_step .AND. (error_max < eps_diis)) THEN
1218 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1219 0 : "DIIS | => A SCF DIIS step will be performed"
1220 : ELSE
1221 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1222 0 : "DIIS | => No SCF DIIS step will be performed"
1223 : END IF
1224 : END IF
1225 :
1226 : ! Update the SCF DIIS buffer
1227 27498 : IF (error_max < eps_diis) THEN
1228 82412 : DO jb = 1, nb
1229 64354 : diis_buffer%b_matrix(ib, jb) = b(ib, jb)
1230 82412 : diis_buffer%b_matrix(jb, ib) = b(jb, ib)
1231 : END DO
1232 : ELSE
1233 :
1234 9440 : diis_step = .FALSE.
1235 : END IF
1236 27498 : DEALLOCATE (b)
1237 :
1238 : ! Perform DIIS step
1239 27498 : IF (diis_step) THEN
1240 :
1241 10694 : nb1 = nb + 1
1242 :
1243 42776 : ALLOCATE (a(nb1, nb1))
1244 32082 : ALLOCATE (b(nb1, nb1))
1245 32082 : ALLOCATE (ev(nb1))
1246 :
1247 : ! Set up the linear DIIS equation system
1248 213814 : b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
1249 :
1250 52004 : b(1:nb, nb1) = -1.0_dp
1251 52004 : b(nb1, 1:nb) = -1.0_dp
1252 10694 : b(nb1, nb1) = 0.0_dp
1253 :
1254 : ! Solve the linear DIIS equation system
1255 62698 : ev(1:nb1) = 0.0_dp !eigenvalues
1256 317822 : a(1:nb1, 1:nb1) = 0.0_dp !eigenvectors
1257 10694 : CALL diag_complex(b(1:nb1, 1:nb1), a(1:nb1, 1:nb1), ev(1:nb1))
1258 317822 : b(1:nb1, 1:nb1) = a(1:nb1, 1:nb1)
1259 :
1260 10694 : eigenvectors_discarded = .FALSE.
1261 :
1262 62698 : DO jb = 1, nb1
1263 62698 : IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
1264 18164 : IF (output_unit > 0) THEN
1265 0 : IF (.NOT. eigenvectors_discarded) THEN
1266 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
1267 0 : "DIIS | Checking eigenvalues of the DIIS error matrix"
1268 : END IF
1269 : WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
1270 0 : "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
1271 0 : "threshold ", eigenvalue_threshold
1272 0 : CALL compress(message)
1273 0 : WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
1274 0 : eigenvectors_discarded = .TRUE.
1275 : END IF
1276 108700 : a(1:nb1, jb) = 0.0_dp
1277 : ELSE
1278 198428 : a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
1279 : END IF
1280 : END DO
1281 :
1282 10694 : IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
1283 : WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
1284 0 : "DIIS | The corresponding eigenvectors were discarded"
1285 : END IF
1286 :
1287 328516 : coeffs(1:nb) = -MATMUL(a(1:nb, 1:nb1), CONJG(b(nb1, 1:nb1)))
1288 : ELSE
1289 :
1290 70790 : coeffs(:) = 0.0_dp
1291 16804 : coeffs(ib) = 1.0_dp
1292 : END IF
1293 :
1294 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1295 27498 : "PRINT%DIIS_INFO")
1296 :
1297 27498 : CALL timestop(handle)
1298 :
1299 82494 : END SUBROUTINE qs_diis_b_step_kp
1300 10694 : END MODULE qs_diis
|
