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