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