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