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 History management and matrix combination for ADIIS.
10 : ! **************************************************************************************************
11 : MODULE qs_scf_subspace
12 :
13 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
14 : dbcsr_copy,&
15 : dbcsr_create,&
16 : dbcsr_p_type,&
17 : dbcsr_set
18 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
19 : USE ieee_arithmetic, ONLY: ieee_is_finite
20 : USE kinds, ONLY: dp
21 : USE qs_scf_subspace_math, ONLY: qs_scf_subspace_build_adiis_model,&
22 : qs_scf_subspace_fifo_slot,&
23 : simplex_qp_nonfinite_input,&
24 : simplex_qp_pairwise_stationary,&
25 : simplex_qp_success,&
26 : simplex_quadratic_minimize
27 : USE qs_scf_subspace_types, ONLY: qs_scf_subspace_buffer_clear,&
28 : qs_scf_subspace_buffer_type
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_subspace'
36 :
37 : INTEGER, PARAMETER, PUBLIC :: scf_subspace_invalid_history = 100
38 :
39 : PUBLIC :: qs_scf_subspace_build, qs_scf_subspace_push, qs_scf_subspace_restart
40 : PUBLIC :: qs_scf_subspace_update_shift
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief Adapt candidate regularization from two evaluated states; no trial KS builds.
46 : !> \param buffer Paired physical states, including the newest evaluated density.
47 : !> \param shift Current occupied-space level shift in hartree.
48 : ! **************************************************************************************************
49 50 : SUBROUTINE qs_scf_subspace_update_shift(buffer, shift)
50 : TYPE(qs_scf_subspace_buffer_type), INTENT(IN) :: buffer
51 : REAL(KIND=dp), INTENT(INOUT) :: shift
52 :
53 : INTEGER :: newest, previous
54 : REAL(KIND=dp) :: actual, predicted, tolerance
55 :
56 : ! Only pure ADIIS steps measure the response to this regularization.
57 50 : IF (buffer%nstored < 2 .OR. buffer%diis_weight > 0.0_dp) RETURN
58 340 : newest = MAXLOC(buffer%generation, DIM=1)
59 : previous = MAXLOC(buffer%generation, DIM=1, &
60 340 : MASK=buffer%generation < buffer%generation(newest))
61 : ! T_ij = Tr(P_i F_j). Compare actual progress with the linear model at P_previous.
62 20 : actual = buffer%state_energy(newest) - buffer%state_energy(previous)
63 20 : predicted = buffer%pf_metric(newest, previous) - buffer%pf_metric(previous, previous)
64 : tolerance = 64.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, &
65 : ABS(buffer%state_energy(newest)), ABS(buffer%state_energy(previous)), &
66 20 : ABS(buffer%pf_metric(newest, previous)), ABS(buffer%pf_metric(previous, previous)))
67 20 : IF (actual > tolerance .OR. predicted > tolerance) THEN
68 : ! An uphill endpoint still contains useful curvature information: retain the history.
69 0 : shift = 2.0_dp*shift
70 20 : ELSE IF (predicted < -tolerance) THEN
71 20 : IF (actual > 0.25_dp*predicted) shift = 2.0_dp*shift
72 20 : IF (actual < 0.75_dp*predicted) shift = 0.5_dp*shift
73 : ELSE
74 : ! Do not maintain a large penalty when changes fall below energy resolution.
75 0 : shift = 0.5_dp*shift
76 : END IF
77 20 : CPASSERT(ieee_is_finite(shift))
78 :
79 : END SUBROUTINE qs_scf_subspace_update_shift
80 :
81 : ! **************************************************************************************************
82 : !> \brief Restart the ADIIS history from the current strictly paired P, F[P] state.
83 : !> \param buffer Persistent SCF subspace history.
84 : !> \param fock Current raw F[P], indexed by spin and real-space cell.
85 : !> \param density Current input P that generated fock, with matching indices.
86 : !> \param state_energy ...
87 : !> \param restarted Whether the current state was stored successfully.
88 : !> \param density_aux Optional replicated density coordinates for a separately represented term.
89 : !> \param fock_aux Conjugate potential for the same state; contributes SUM(density_aux*fock_aux) to Tr(PF).
90 : ! **************************************************************************************************
91 0 : SUBROUTINE qs_scf_subspace_restart(buffer, fock, density, state_energy, restarted, density_aux, fock_aux)
92 :
93 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
94 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: fock, density
95 : REAL(KIND=dp), INTENT(IN) :: state_energy
96 : LOGICAL, INTENT(OUT) :: restarted
97 : REAL(dp), INTENT(IN), OPTIONAL :: density_aux(:, :), fock_aux(:, :)
98 :
99 0 : restarted = .FALSE.
100 0 : CALL qs_scf_subspace_buffer_clear(buffer)
101 0 : CALL qs_scf_subspace_push(buffer, fock, density, state_energy, restarted, density_aux, fock_aux)
102 0 : IF (.NOT. restarted) RETURN
103 :
104 0 : buffer%coefficients = 0.0_dp
105 0 : buffer%coefficients(1) = 1.0_dp
106 0 : buffer%last_status = simplex_qp_success
107 0 : buffer%last_restart = .TRUE.
108 0 : buffer%use_combined_fock = .FALSE.
109 0 : buffer%last_objective = 0.0_dp
110 0 : buffer%last_old_fock_weight = 0.0_dp
111 :
112 : END SUBROUTINE qs_scf_subspace_restart
113 :
114 : ! **************************************************************************************************
115 : !> \brief Append one accepted, evaluated P,F[P] state to the ADIIS history.
116 : !>
117 : !> The SCF driver owns state acceptance. Trial endpoints and arbitrary initial guesses
118 : !> must not call this routine.
119 : !> \param buffer Persistent SCF subspace history.
120 : !> \param fock Current raw F[P], indexed by spin and real-space cell.
121 : !> \param density Current input P that generated fock, with matching indices.
122 : !> \param state_energy ...
123 : !> \param pushed Whether the state was stored successfully.
124 : !> \param density_aux Optional replicated density coordinates for a separately represented term.
125 : !> \param fock_aux Conjugate potential for the same state; contributes SUM(density_aux*fock_aux) to Tr(PF).
126 : ! **************************************************************************************************
127 110 : SUBROUTINE qs_scf_subspace_push(buffer, fock, density, state_energy, pushed, density_aux, fock_aux)
128 :
129 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
130 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: fock, density
131 : REAL(KIND=dp), INTENT(IN) :: state_energy
132 : LOGICAL, INTENT(OUT) :: pushed
133 : REAL(dp), INTENT(IN), OPTIONAL :: density_aux(:, :), fock_aux(:, :)
134 :
135 : LOGICAL :: storage_valid
136 :
137 110 : pushed = .FALSE.
138 110 : buffer%use_combined_fock = .FALSE.
139 110 : buffer%last_restart = .FALSE.
140 110 : buffer%last_objective = 0.0_dp
141 110 : buffer%last_old_fock_weight = 0.0_dp
142 1678 : IF (ALLOCATED(buffer%coefficients)) buffer%coefficients = 0.0_dp
143 :
144 110 : IF (buffer%nbuffer < 1) THEN
145 0 : buffer%last_status = scf_subspace_invalid_history
146 0 : RETURN
147 : END IF
148 110 : IF (buffer%ncall < 0 .OR. buffer%nstored < 0 .OR. buffer%nstored > buffer%nbuffer .OR. &
149 : buffer%nstored /= MIN(buffer%ncall, buffer%nbuffer)) THEN
150 0 : buffer%last_status = scf_subspace_invalid_history
151 0 : RETURN
152 : END IF
153 :
154 110 : CALL ensure_matrix_storage(buffer, fock, density, storage_valid)
155 110 : IF (.NOT. storage_valid) THEN
156 0 : buffer%last_status = scf_subspace_invalid_history
157 0 : RETURN
158 : END IF
159 :
160 110 : CPASSERT(PRESENT(density_aux) .EQV. PRESENT(fock_aux))
161 110 : IF (PRESENT(density_aux)) THEN
162 0 : CPASSERT(ALL(SHAPE(density_aux) == SHAPE(fock_aux)))
163 0 : IF (.NOT. ALLOCATED(buffer%density_aux)) THEN
164 0 : CPASSERT(buffer%nstored == 0)
165 0 : ALLOCATE (buffer%density_aux(SIZE(density_aux, 1), SIZE(density_aux, 2), buffer%nbuffer), &
166 0 : buffer%fock_aux(SIZE(density_aux, 1), SIZE(density_aux, 2), buffer%nbuffer), &
167 0 : buffer%combined_aux(SIZE(density_aux, 1), SIZE(density_aux, 2)))
168 : END IF
169 0 : CPASSERT(ALL(SHAPE(buffer%combined_aux) == SHAPE(density_aux)))
170 : ELSE
171 110 : CPASSERT(.NOT. ALLOCATED(buffer%density_aux))
172 : END IF
173 330 : CALL push_paired_state(buffer, fock, density, state_energy, pushed, density_aux, fock_aux)
174 110 : IF (.NOT. pushed) THEN
175 0 : buffer%last_status = simplex_qp_nonfinite_input
176 0 : RETURN
177 : END IF
178 110 : buffer%last_status = simplex_qp_success
179 :
180 : END SUBROUTINE qs_scf_subspace_push
181 :
182 : ! **************************************************************************************************
183 : !> \brief Solve the ADIIS model from previously accepted history and form an effective KS matrix.
184 : !> \param buffer Persistent SCF subspace history.
185 : ! **************************************************************************************************
186 110 : SUBROUTINE qs_scf_subspace_build(buffer)
187 :
188 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
189 :
190 : INTEGER :: i, j, m, physical
191 110 : INTEGER, ALLOCATABLE, DIMENSION(:) :: slots
192 : LOGICAL :: model_valid, solver_usable
193 : REAL(KIND=dp) :: coefficient_sum, model_delta, &
194 : model_tolerance
195 110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: difference, linear
196 110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chronological_pf, hessian
197 :
198 110 : buffer%use_combined_fock = .FALSE.
199 110 : buffer%last_restart = .FALSE.
200 110 : buffer%last_objective = 0.0_dp
201 110 : buffer%last_old_fock_weight = 0.0_dp
202 1870 : IF (ALLOCATED(buffer%coefficients)) buffer%coefficients = 0.0_dp
203 :
204 110 : IF (buffer%nbuffer < 1 .OR. buffer%nstored < 1) THEN
205 0 : buffer%last_status = scf_subspace_invalid_history
206 0 : RETURN
207 : END IF
208 110 : IF (buffer%ncall < 0 .OR. buffer%nstored < 0 .OR. buffer%nstored > buffer%nbuffer .OR. &
209 : buffer%nstored /= MIN(buffer%ncall, buffer%nbuffer)) THEN
210 0 : buffer%last_status = scf_subspace_invalid_history
211 0 : RETURN
212 : END IF
213 :
214 110 : IF (.NOT. ASSOCIATED(buffer%fock) .OR. .NOT. ASSOCIATED(buffer%density) .OR. &
215 : .NOT. ASSOCIATED(buffer%combined_fock)) THEN
216 0 : buffer%last_status = scf_subspace_invalid_history
217 0 : RETURN
218 : END IF
219 : IF (.NOT. ALLOCATED(buffer%generation) .OR. .NOT. ALLOCATED(buffer%pf_metric) .OR. &
220 110 : .NOT. ALLOCATED(buffer%coefficients) .OR. .NOT. ALLOCATED(buffer%state_energy)) THEN
221 0 : buffer%last_status = scf_subspace_invalid_history
222 0 : RETURN
223 : END IF
224 : IF (SIZE(buffer%generation) /= buffer%nbuffer .OR. &
225 : SIZE(buffer%pf_metric, 1) /= buffer%nbuffer .OR. SIZE(buffer%pf_metric, 2) /= buffer%nbuffer .OR. &
226 : SIZE(buffer%coefficients) /= buffer%nbuffer .OR. SIZE(buffer%state_energy) /= buffer%nbuffer .OR. &
227 3630 : COUNT(buffer%generation > 0) /= buffer%nstored .OR. ANY(buffer%generation < 0)) THEN
228 0 : buffer%last_status = scf_subspace_invalid_history
229 0 : RETURN
230 : END IF
231 :
232 110 : m = buffer%nstored
233 1100 : ALLOCATE (slots(m), chronological_pf(m, m), hessian(m, m), linear(m))
234 110 : CALL chronological_slots(buffer, slots)
235 :
236 720 : DO i = 1, m
237 5266 : DO j = 1, m
238 5156 : chronological_pf(i, j) = buffer%pf_metric(slots(i), slots(j))
239 : END DO
240 : END DO
241 :
242 110 : CALL qs_scf_subspace_build_adiis_model(chronological_pf, m, hessian, linear, model_valid)
243 :
244 1870 : buffer%coefficients = 0.0_dp
245 110 : IF (.NOT. model_valid) THEN
246 0 : CALL restart_newest_after_failure(buffer, scf_subspace_invalid_history)
247 0 : RETURN
248 : END IF
249 :
250 : CALL simplex_quadratic_minimize(hessian, linear, buffer%coefficients(1:m), &
251 110 : buffer%last_objective, buffer%last_status, preferred_index=m)
252 :
253 : solver_usable = buffer%last_status == simplex_qp_success .OR. &
254 110 : buffer%last_status == simplex_qp_pairwise_stationary
255 110 : IF (.NOT. solver_usable) THEN
256 0 : CALL restart_newest_after_failure(buffer, buffer%last_status)
257 0 : RETURN
258 : END IF
259 720 : IF (.NOT. ieee_is_finite(buffer%last_objective) .OR. &
260 : .NOT. ALL(ieee_is_finite(buffer%coefficients(1:m)))) THEN
261 0 : CALL restart_newest_after_failure(buffer, simplex_qp_nonfinite_input)
262 0 : RETURN
263 : END IF
264 :
265 720 : coefficient_sum = SUM(buffer%coefficients(1:m))
266 720 : IF (MINVAL(buffer%coefficients(1:m)) < -100.0_dp*EPSILON(1.0_dp) .OR. &
267 : ABS(coefficient_sum - 1.0_dp) > 1000.0_dp*EPSILON(1.0_dp)) THEN
268 0 : CALL restart_newest_after_failure(buffer, scf_subspace_invalid_history)
269 0 : RETURN
270 : END IF
271 :
272 220 : ALLOCATE (difference(m))
273 720 : difference(:) = buffer%coefficients(1:m)
274 110 : difference(m) = difference(m) - 1.0_dp
275 330 : model_delta = 0.5_dp*DOT_PRODUCT(difference, MATMUL(hessian, difference)) + &
276 11032 : DOT_PRODUCT(difference, hessian(:, m) + linear)
277 : model_tolerance = 1000.0_dp*EPSILON(1.0_dp)*REAL(m, KIND=dp)* &
278 5876 : MAX(1.0_dp, MAXVAL(ABS(hessian)), MAXVAL(ABS(linear)))
279 110 : IF (.NOT. ieee_is_finite(model_delta) .OR. model_delta > model_tolerance) THEN
280 0 : IF (ieee_is_finite(model_delta)) THEN
281 0 : CALL restart_newest_after_failure(buffer, scf_subspace_invalid_history)
282 : ELSE
283 0 : CALL restart_newest_after_failure(buffer, simplex_qp_nonfinite_input)
284 : END IF
285 0 : RETURN
286 : END IF
287 :
288 110 : buffer%last_old_fock_weight = MAX(0.0_dp, 1.0_dp - buffer%coefficients(m))
289 110 : IF (buffer%last_old_fock_weight <= 100.0_dp*EPSILON(1.0_dp)) RETURN
290 :
291 62 : CALL zero_fock_combination(buffer)
292 340 : DO j = 1, m
293 278 : physical = slots(j)
294 278 : IF (buffer%coefficients(j) <= 0.0_dp) CYCLE
295 340 : CALL add_fock_to_combination(buffer, physical, buffer%coefficients(j))
296 : END DO
297 62 : buffer%use_combined_fock = .TRUE.
298 :
299 220 : END SUBROUTINE qs_scf_subspace_build
300 :
301 : ! **************************************************************************************************
302 : !> \brief Retain only the newest accepted state after an unusable ADIIS model.
303 : !> \param buffer Persistent SCF subspace history.
304 : !> \param failure_status Status code that caused the restart.
305 : ! **************************************************************************************************
306 0 : SUBROUTINE restart_newest_after_failure(buffer, failure_status)
307 :
308 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
309 : INTEGER, INTENT(IN) :: failure_status
310 :
311 : INTEGER :: icell, ispin, newest
312 : REAL(KIND=dp) :: contribution, newest_energy, newest_pf
313 :
314 0 : IF (buffer%nstored < 1) THEN
315 0 : buffer%last_status = failure_status
316 0 : RETURN
317 : END IF
318 0 : IF (COUNT(buffer%generation > 0) /= buffer%nstored) THEN
319 0 : buffer%last_status = failure_status
320 0 : RETURN
321 : END IF
322 0 : newest = MAXLOC(buffer%generation, DIM=1)
323 0 : newest_energy = buffer%state_energy(newest)
324 0 : newest_pf = 0.0_dp
325 0 : DO icell = 1, SIZE(buffer%fock, 3)
326 0 : DO ispin = 1, SIZE(buffer%fock, 2)
327 : CALL dbcsr_dot(buffer%density(newest, ispin, icell)%matrix, &
328 0 : buffer%fock(newest, ispin, icell)%matrix, contribution)
329 0 : newest_pf = newest_pf + contribution
330 : END DO
331 : END DO
332 0 : IF (ALLOCATED(buffer%density_aux)) THEN
333 0 : newest_pf = newest_pf + SUM(buffer%density_aux(:, :, newest)*buffer%fock_aux(:, :, newest))
334 : END IF
335 0 : IF (.NOT. ieee_is_finite(newest_energy) .OR. .NOT. ieee_is_finite(newest_pf)) THEN
336 0 : buffer%last_status = failure_status
337 0 : RETURN
338 : END IF
339 :
340 0 : buffer%ncall = 1
341 0 : buffer%nstored = 1
342 0 : buffer%generation = 0
343 0 : buffer%generation(newest) = 1
344 0 : buffer%state_energy = HUGE(1.0_dp)
345 0 : buffer%state_energy(newest) = newest_energy
346 0 : buffer%pf_metric = 0.0_dp
347 0 : buffer%pf_metric(newest, newest) = newest_pf
348 0 : buffer%coefficients = 0.0_dp
349 0 : buffer%coefficients(1) = 1.0_dp
350 0 : buffer%last_status = failure_status
351 0 : buffer%last_restart = .TRUE.
352 0 : buffer%use_combined_fock = .FALSE.
353 0 : buffer%last_objective = 0.0_dp
354 0 : buffer%last_old_fock_weight = 0.0_dp
355 0 : buffer%diis_weight = 0.0_dp
356 :
357 : END SUBROUTINE restart_newest_after_failure
358 :
359 : ! **************************************************************************************************
360 : !> \brief Allocate matrix and scalar storage lazily from current matrix templates.
361 : !> \param buffer Persistent SCF subspace history.
362 : !> \param fock Current Fock-matrix templates.
363 : !> \param density Current density-matrix templates.
364 : !> \param valid Whether dimensions and existing allocation match.
365 : ! **************************************************************************************************
366 110 : SUBROUTINE ensure_matrix_storage(buffer, fock, density, valid)
367 :
368 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
369 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: fock, density
370 : LOGICAL, INTENT(OUT) :: valid
371 :
372 : INTEGER :: icell, islot, ispin, ncell, nspin
373 :
374 110 : valid = ASSOCIATED(fock) .AND. ASSOCIATED(density)
375 110 : IF (.NOT. valid) RETURN
376 110 : nspin = SIZE(fock, 1)
377 110 : ncell = SIZE(fock, 2)
378 : valid = nspin > 0 .AND. ncell > 0 .AND. &
379 110 : SIZE(density, 1) == nspin .AND. SIZE(density, 2) == ncell
380 110 : IF (.NOT. valid) RETURN
381 6042 : DO icell = 1, ncell
382 17258 : DO ispin = 1, nspin
383 : valid = ASSOCIATED(fock(ispin, icell)%matrix) .AND. &
384 11216 : ASSOCIATED(density(ispin, icell)%matrix)
385 17148 : IF (.NOT. valid) RETURN
386 : END DO
387 : END DO
388 :
389 110 : IF (ASSOCIATED(buffer%fock)) THEN
390 : valid = SIZE(buffer%fock, 1) == buffer%nbuffer .AND. &
391 98 : SIZE(buffer%fock, 2) == nspin .AND. SIZE(buffer%fock, 3) == ncell
392 98 : IF (.NOT. valid) RETURN
393 98 : valid = ASSOCIATED(buffer%density)
394 98 : IF (.NOT. valid) RETURN
395 : valid = SIZE(buffer%density, 1) == buffer%nbuffer .AND. &
396 98 : SIZE(buffer%density, 2) == nspin .AND. SIZE(buffer%density, 3) == ncell
397 98 : IF (.NOT. valid) RETURN
398 98 : valid = ASSOCIATED(buffer%combined_fock)
399 98 : IF (.NOT. valid) RETURN
400 : valid = SIZE(buffer%combined_fock, 1) == nspin .AND. &
401 98 : SIZE(buffer%combined_fock, 2) == ncell
402 98 : IF (.NOT. valid) RETURN
403 : valid = ALLOCATED(buffer%generation) .AND. ALLOCATED(buffer%pf_metric) .AND. &
404 98 : ALLOCATED(buffer%coefficients) .AND. ALLOCATED(buffer%state_energy)
405 98 : IF (.NOT. valid) RETURN
406 : valid = SIZE(buffer%generation) == buffer%nbuffer .AND. &
407 : SIZE(buffer%pf_metric, 1) == buffer%nbuffer .AND. &
408 : SIZE(buffer%pf_metric, 2) == buffer%nbuffer .AND. &
409 : SIZE(buffer%coefficients) == buffer%nbuffer .AND. &
410 98 : SIZE(buffer%state_energy) == buffer%nbuffer
411 98 : IF (.NOT. valid) RETURN
412 5476 : DO icell = 1, ncell
413 15656 : DO ispin = 1, nspin
414 10180 : valid = ASSOCIATED(buffer%combined_fock(ispin, icell)%matrix)
415 10180 : IF (.NOT. valid) RETURN
416 178438 : DO islot = 1, buffer%nbuffer
417 : valid = ASSOCIATED(buffer%fock(islot, ispin, icell)%matrix) .AND. &
418 162880 : ASSOCIATED(buffer%density(islot, ispin, icell)%matrix)
419 173060 : IF (.NOT. valid) RETURN
420 : END DO
421 : END DO
422 : END DO
423 : RETURN
424 : END IF
425 :
426 : valid = .NOT. ASSOCIATED(buffer%density) .AND. &
427 : .NOT. ASSOCIATED(buffer%combined_fock) .AND. .NOT. ALLOCATED(buffer%generation) .AND. &
428 : .NOT. ALLOCATED(buffer%pf_metric) .AND. .NOT. ALLOCATED(buffer%coefficients) .AND. &
429 12 : .NOT. ALLOCATED(buffer%state_energy)
430 12 : IF (.NOT. valid) RETURN
431 :
432 18226 : ALLOCATE (buffer%fock(buffer%nbuffer, nspin, ncell))
433 18214 : ALLOCATE (buffer%density(buffer%nbuffer, nspin, ncell))
434 1638 : ALLOCATE (buffer%combined_fock(nspin, ncell))
435 60 : ALLOCATE (buffer%generation(buffer%nbuffer), buffer%state_energy(buffer%nbuffer))
436 60 : ALLOCATE (buffer%pf_metric(buffer%nbuffer, buffer%nbuffer), buffer%coefficients(buffer%nbuffer))
437 204 : buffer%generation = 0
438 204 : buffer%state_energy = HUGE(1.0_dp)
439 3276 : buffer%pf_metric = 0.0_dp
440 204 : buffer%coefficients = 0.0_dp
441 :
442 566 : DO icell = 1, ncell
443 1602 : DO ispin = 1, nspin
444 17612 : DO islot = 1, buffer%nbuffer
445 16576 : ALLOCATE (buffer%fock(islot, ispin, icell)%matrix)
446 : CALL dbcsr_create(buffer%fock(islot, ispin, icell)%matrix, &
447 16576 : template=fock(ispin, icell)%matrix)
448 16576 : ALLOCATE (buffer%density(islot, ispin, icell)%matrix)
449 : CALL dbcsr_create(buffer%density(islot, ispin, icell)%matrix, &
450 17612 : template=density(ispin, icell)%matrix)
451 : END DO
452 1036 : ALLOCATE (buffer%combined_fock(ispin, icell)%matrix)
453 : CALL dbcsr_create(buffer%combined_fock(ispin, icell)%matrix, &
454 1590 : template=fock(ispin, icell)%matrix)
455 : END DO
456 : END DO
457 :
458 : END SUBROUTINE ensure_matrix_storage
459 :
460 : ! **************************************************************************************************
461 : !> \brief Atomically append one strictly paired P, F[P] state.
462 : !> \param buffer Persistent SCF subspace history.
463 : !> \param fock Current raw F[P].
464 : !> \param density Current P.
465 : !> \param state_energy ...
466 : !> \param pushed Whether all scalar checks passed and the state was committed.
467 : !> \param density_aux Optional replicated density coordinates for a separately represented term.
468 : !> \param fock_aux Conjugate potential for the same state; contributes SUM(density_aux*fock_aux) to Tr(PF).
469 : ! **************************************************************************************************
470 110 : SUBROUTINE push_paired_state(buffer, fock, density, state_energy, pushed, density_aux, fock_aux)
471 :
472 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
473 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: fock, density
474 : REAL(KIND=dp), INTENT(IN) :: state_energy
475 : LOGICAL, INTENT(OUT) :: pushed
476 : REAL(dp), INTENT(IN), OPTIONAL :: density_aux(:, :), fock_aux(:, :)
477 :
478 : INTEGER :: icell, islot, ispin, target_slot
479 : LOGICAL :: retained
480 : REAL(KIND=dp) :: current_current
481 110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: current_old, old_current
482 :
483 110 : pushed = .FALSE.
484 110 : IF (.NOT. ieee_is_finite(state_energy)) RETURN
485 :
486 : ! Keep the baseline deterministic and auditable: fill empty slots in physical
487 : ! order, then evict the oldest accepted state. Energy remains metadata.
488 110 : target_slot = qs_scf_subspace_fifo_slot(buffer%generation)
489 110 : IF (target_slot <= 0) RETURN
490 :
491 440 : ALLOCATE (current_old(buffer%nbuffer), old_current(buffer%nbuffer))
492 110 : current_old = 0.0_dp
493 110 : old_current = 0.0_dp
494 :
495 110 : CALL paired_matrix_dot(density, fock, current_current)
496 110 : IF (PRESENT(density_aux)) current_current = current_current + SUM(density_aux*fock_aux)
497 110 : IF (.NOT. ieee_is_finite(current_current)) RETURN
498 :
499 1870 : DO islot = 1, buffer%nbuffer
500 2260 : retained = buffer%generation(islot) > 0 .AND. islot /= target_slot
501 : IF (.NOT. retained) CYCLE
502 500 : CALL history_current_dot(buffer, islot, fock, density, old_current(islot), current_old(islot))
503 500 : IF (PRESENT(density_aux)) THEN
504 0 : old_current(islot) = old_current(islot) + SUM(buffer%density_aux(:, :, islot)*fock_aux)
505 0 : current_old(islot) = current_old(islot) + SUM(density_aux*buffer%fock_aux(:, :, islot))
506 : END IF
507 500 : IF (.NOT. ieee_is_finite(old_current(islot)) .OR. &
508 110 : .NOT. ieee_is_finite(current_old(islot))) RETURN
509 : END DO
510 :
511 6042 : DO icell = 1, SIZE(fock, 2)
512 17258 : DO ispin = 1, SIZE(fock, 1)
513 11216 : CALL dbcsr_copy(buffer%fock(target_slot, ispin, icell)%matrix, fock(ispin, icell)%matrix)
514 17148 : CALL dbcsr_copy(buffer%density(target_slot, ispin, icell)%matrix, density(ispin, icell)%matrix)
515 : END DO
516 : END DO
517 :
518 110 : IF (PRESENT(density_aux)) THEN
519 0 : buffer%density_aux(:, :, target_slot) = density_aux
520 0 : buffer%fock_aux(:, :, target_slot) = fock_aux
521 : END IF
522 1870 : buffer%pf_metric(target_slot, :) = 0.0_dp
523 1870 : buffer%pf_metric(:, target_slot) = 0.0_dp
524 110 : buffer%pf_metric(target_slot, target_slot) = current_current
525 1870 : DO islot = 1, buffer%nbuffer
526 2260 : retained = buffer%generation(islot) > 0 .AND. islot /= target_slot
527 : IF (.NOT. retained) CYCLE
528 500 : buffer%pf_metric(islot, target_slot) = old_current(islot)
529 610 : buffer%pf_metric(target_slot, islot) = current_old(islot)
530 : END DO
531 :
532 110 : buffer%ncall = buffer%ncall + 1
533 110 : buffer%generation(target_slot) = buffer%ncall
534 110 : buffer%state_energy(target_slot) = state_energy
535 110 : buffer%nstored = MIN(buffer%nstored + 1, buffer%nbuffer)
536 110 : pushed = .TRUE.
537 :
538 220 : END SUBROUTINE push_paired_state
539 :
540 : ! **************************************************************************************************
541 : !> \brief Compute Tr(P F), summed across every spin and real-space cell.
542 : !> \param density Density matrices.
543 : !> \param fock Fock matrices.
544 : !> \param value Summed trace contraction.
545 : ! **************************************************************************************************
546 110 : SUBROUTINE paired_matrix_dot(density, fock, value)
547 :
548 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density, fock
549 : REAL(KIND=dp), INTENT(OUT) :: value
550 :
551 : INTEGER :: icell, ispin
552 : REAL(KIND=dp) :: contribution
553 :
554 110 : value = 0.0_dp
555 6042 : DO icell = 1, SIZE(fock, 2)
556 17258 : DO ispin = 1, SIZE(fock, 1)
557 11216 : CALL dbcsr_dot(density(ispin, icell)%matrix, fock(ispin, icell)%matrix, contribution)
558 17148 : value = value + contribution
559 : END DO
560 : END DO
561 :
562 110 : END SUBROUTINE paired_matrix_dot
563 :
564 : ! **************************************************************************************************
565 : !> \brief Compute both cross contractions between one old and the current state.
566 : !> \param buffer Persistent SCF subspace history.
567 : !> \param islot Physical old-history slot.
568 : !> \param current_fock Current F[P].
569 : !> \param current_density Current P.
570 : !> \param old_current Tr(P_old F_current).
571 : !> \param current_old Tr(P_current F_old).
572 : ! **************************************************************************************************
573 500 : SUBROUTINE history_current_dot(buffer, islot, current_fock, current_density, old_current, current_old)
574 :
575 : TYPE(qs_scf_subspace_buffer_type), INTENT(IN) :: buffer
576 : INTEGER, INTENT(IN) :: islot
577 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: current_fock, current_density
578 : REAL(KIND=dp), INTENT(OUT) :: old_current, current_old
579 :
580 : INTEGER :: icell, ispin
581 : REAL(KIND=dp) :: contribution
582 :
583 500 : old_current = 0.0_dp
584 500 : current_old = 0.0_dp
585 30424 : DO icell = 1, SIZE(current_fock, 2)
586 87680 : DO ispin = 1, SIZE(current_fock, 1)
587 : CALL dbcsr_dot(buffer%density(islot, ispin, icell)%matrix, &
588 57256 : current_fock(ispin, icell)%matrix, contribution)
589 57256 : old_current = old_current + contribution
590 : CALL dbcsr_dot(current_density(ispin, icell)%matrix, &
591 57256 : buffer%fock(islot, ispin, icell)%matrix, contribution)
592 87180 : current_old = current_old + contribution
593 : END DO
594 : END DO
595 :
596 500 : END SUBROUTINE history_current_dot
597 :
598 : ! **************************************************************************************************
599 : !> \brief Return physical slots ordered from oldest to newest.
600 : !> \param buffer Persistent SCF subspace history.
601 : !> \param slots Chronologically ordered physical slot indices.
602 : ! **************************************************************************************************
603 110 : PURE SUBROUTINE chronological_slots(buffer, slots)
604 :
605 : TYPE(qs_scf_subspace_buffer_type), INTENT(IN) :: buffer
606 : INTEGER, DIMENSION(:), INTENT(OUT) :: slots
607 :
608 : INTEGER :: i, islot, j, tmp
609 :
610 110 : j = 0
611 1870 : DO islot = 1, buffer%nbuffer
612 1760 : IF (buffer%generation(islot) <= 0) CYCLE
613 610 : j = j + 1
614 1870 : IF (j <= SIZE(slots)) slots(j) = islot
615 : END DO
616 :
617 : ! The ADIIS model uses the most recent state as its expansion point, so
618 : ! return physical slots ordered by insertion generation with the newest
619 : ! state last. The history size is small, making insertion sort sufficient.
620 610 : DO i = 2, SIZE(slots)
621 500 : tmp = slots(i)
622 500 : j = i - 1
623 500 : DO WHILE (j >= 1)
624 500 : IF (buffer%generation(slots(j)) <= buffer%generation(tmp)) EXIT
625 0 : slots(j + 1) = slots(j)
626 500 : j = j - 1
627 : END DO
628 610 : slots(j + 1) = tmp
629 : END DO
630 :
631 110 : END SUBROUTINE chronological_slots
632 :
633 : ! **************************************************************************************************
634 : !> \brief Zero the effective KS-matrix combination.
635 : !> \param buffer Persistent SCF subspace history and combination target.
636 : ! **************************************************************************************************
637 62 : SUBROUTINE zero_fock_combination(buffer)
638 :
639 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
640 :
641 : INTEGER :: icell, ispin
642 :
643 62 : IF (ALLOCATED(buffer%combined_aux)) buffer%combined_aux(:, :) = 0.0_dp
644 2904 : DO icell = 1, SIZE(buffer%combined_fock, 2)
645 8156 : DO ispin = 1, SIZE(buffer%combined_fock, 1)
646 8094 : CALL dbcsr_set(buffer%combined_fock(ispin, icell)%matrix, 0.0_dp)
647 : END DO
648 : END DO
649 :
650 62 : END SUBROUTINE zero_fock_combination
651 :
652 : ! **************************************************************************************************
653 : !> \brief Add one physical history slot to the effective KS-matrix combination.
654 : !> \param buffer Persistent SCF subspace history and combination target.
655 : !> \param physical Physical history slot.
656 : !> \param coefficient Weight of this history matrix.
657 : ! **************************************************************************************************
658 124 : SUBROUTINE add_fock_to_combination(buffer, physical, coefficient)
659 :
660 : TYPE(qs_scf_subspace_buffer_type), INTENT(INOUT) :: buffer
661 : INTEGER, INTENT(IN) :: physical
662 : REAL(KIND=dp), INTENT(IN) :: coefficient
663 :
664 : INTEGER :: icell, ispin
665 :
666 124 : IF (ALLOCATED(buffer%combined_aux)) THEN
667 0 : buffer%combined_aux(:, :) = buffer%combined_aux + coefficient*buffer%fock_aux(:, :, physical)
668 : END IF
669 5808 : DO icell = 1, SIZE(buffer%combined_fock, 2)
670 16312 : DO ispin = 1, SIZE(buffer%combined_fock, 1)
671 : CALL dbcsr_add(buffer%combined_fock(ispin, icell)%matrix, &
672 : buffer%fock(physical, ispin, icell)%matrix, &
673 16188 : alpha_scalar=1.0_dp, beta_scalar=coefficient)
674 : END DO
675 : END DO
676 :
677 124 : END SUBROUTINE add_fock_to_combination
678 :
679 110 : END MODULE qs_scf_subspace
|
