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 Safeguarded optimal damping of raw Roothaan SCF steps.
10 : ! **************************************************************************************************
11 : MODULE qs_scf_oda
12 :
13 : USE bibliography, ONLY: Cances2000,&
14 : Herbst2022,&
15 : cite_reference
16 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
17 : dbcsr_copy,&
18 : dbcsr_create,&
19 : dbcsr_get_info,&
20 : dbcsr_p_type
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
22 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
23 : dbcsr_deallocate_matrix_set
24 : USE ieee_arithmetic, ONLY: ieee_is_finite
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE qs_energy_types, ONLY: qs_energy_type
28 : USE qs_environment_types, ONLY: get_qs_env,&
29 : qs_environment_type
30 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
31 : USE qs_ks_types, ONLY: qs_ks_env_type
32 : USE qs_rho_types, ONLY: qs_rho_type
33 : USE qs_scf_loop_utils, ONLY: qs_scf_rho_update
34 : USE qs_scf_types, ONLY: qs_scf_env_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_oda'
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: status_no_acceptable = "NO_ACCEPTABLE", &
43 : status_non_descent = "NON_DESCENT", &
44 : status_nonfinite = "NONFINITE", &
45 : status_success = "SUCCESS"
46 :
47 : PUBLIC :: qs_scf_oda_apply
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief Apply a safeguarded ODA line search to an already diagonalized raw SCF endpoint.
53 : !> \param qs_env Quickstep environment.
54 : !> \param scf_env SCF environment; p_mix_new contains the raw endpoint.
55 : !> \param rho Current accepted density at lambda=0.
56 : !> \param ks_env Kohn-Sham environment.
57 : !> \param base_fock Raw F[P0] retained with the accepted ADIIS history state.
58 : !> \param matrix_ks Fock matrix updated during trial evaluations.
59 : !> \param rho_ao Current accepted AO density, indexed by spin and real-space image cell.
60 : !> \param base_energy Evaluated energy components of P0.
61 : !> \param predictor_lambda Initial trial step and prediction for the next consecutive ODA call.
62 : !> \param predictor_valid Whether predictor_lambda can be reused for a consecutive ODA call.
63 : !> \param applied Whether an energy-safeguarded ODA step was accepted.
64 : !> \param state_evaluated Whether the returned density and Fock state were fully evaluated.
65 : !> \param evaluated_energy Energy components of the fully evaluated returned state.
66 : ! **************************************************************************************************
67 34 : SUBROUTINE qs_scf_oda_apply(qs_env, scf_env, rho, ks_env, base_fock, matrix_ks, rho_ao, base_energy, &
68 : predictor_lambda, predictor_valid, applied, state_evaluated, &
69 : evaluated_energy)
70 :
71 : TYPE(qs_environment_type), POINTER :: qs_env
72 : TYPE(qs_scf_env_type), POINTER :: scf_env
73 : TYPE(qs_rho_type), POINTER :: rho
74 : TYPE(qs_ks_env_type), POINTER :: ks_env
75 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: base_fock
76 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, rho_ao
77 : TYPE(qs_energy_type), INTENT(IN) :: base_energy
78 : REAL(KIND=dp), INTENT(INOUT) :: predictor_lambda
79 : LOGICAL, INTENT(INOUT) :: predictor_valid
80 : LOGICAL, INTENT(OUT) :: applied, state_evaluated
81 : TYPE(qs_energy_type), INTENT(OUT) :: evaluated_energy
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_oda_apply'
84 : INTEGER, PARAMETER :: max_backtracking = 9
85 : REAL(KIND=dp), PARAMETER :: armijo_factor = 1.0E-4_dp, initial_lambda = 0.5_dp, &
86 : maximum_backtrack_fraction = 0.8_dp, minimum_backtrack_fraction = 0.1_dp, &
87 : minimum_lambda = 1.0_dp/1024.0_dp, predictor_growth = 1.5_dp
88 :
89 : CHARACTER(len=15) :: status
90 : CHARACTER(len=default_string_length) :: name
91 : INTEGER :: evaluations, handle, i, icell, ispin
92 : LOGICAL :: model_valid, trial_accepted
93 : REAL(KIND=dp) :: current_lambda, gradient0, gradient1, &
94 : lambda, proposed_fraction, &
95 : proposed_lambda, trial_energy, &
96 : trial_lambda
97 34 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: direction
98 : TYPE(qs_energy_type), POINTER :: energy
99 :
100 34 : CALL timeset(routineN, handle)
101 34 : CALL cite_reference(Cances2000)
102 34 : CALL cite_reference(Herbst2022)
103 34 : NULLIFY (direction, energy)
104 34 : CALL get_qs_env(qs_env, energy=energy)
105 :
106 34 : applied = .FALSE.
107 34 : state_evaluated = .FALSE.
108 34 : lambda = 1.0_dp
109 34 : trial_energy = base_energy%total
110 : gradient0 = 0.0_dp
111 34 : gradient1 = 0.0_dp
112 34 : evaluations = 0
113 34 : status = status_nonfinite
114 :
115 : IF (predictor_valid .AND. ieee_is_finite(predictor_lambda) .AND. &
116 34 : predictor_lambda >= minimum_lambda .AND. predictor_lambda <= 1.0_dp) THEN
117 26 : trial_lambda = predictor_lambda
118 : ELSE
119 8 : trial_lambda = initial_lambda
120 : END IF
121 34 : predictor_lambda = 1.0_dp
122 34 : predictor_valid = .FALSE.
123 :
124 34 : CPASSERT(ASSOCIATED(energy))
125 34 : evaluated_energy = base_energy
126 34 : CPASSERT(ASSOCIATED(scf_env%p_mix_new))
127 34 : CPASSERT(SIZE(scf_env%p_mix_new, 1) == SIZE(rho_ao, 1))
128 34 : CPASSERT(SIZE(scf_env%p_mix_new, 2) == SIZE(rho_ao, 2))
129 34 : CPASSERT(SIZE(base_fock, 1) == SIZE(rho_ao, 1))
130 34 : CPASSERT(SIZE(base_fock, 2) == SIZE(rho_ao, 2))
131 34 : CPASSERT(SIZE(matrix_ks, 1) == SIZE(rho_ao, 1))
132 34 : CPASSERT(SIZE(matrix_ks, 2) == SIZE(rho_ao, 2))
133 :
134 : oda_search: BLOCK
135 34 : CALL create_direction(scf_env%p_mix_new, rho_ao, direction)
136 34 : CALL directional_derivative(base_fock, direction, gradient0)
137 34 : IF (.NOT. ieee_is_finite(gradient0)) EXIT oda_search
138 34 : IF (gradient0 >= 0.0_dp) THEN
139 0 : status = status_non_descent
140 : EXIT oda_search
141 : END IF
142 :
143 46 : DO i = 0, max_backtracking
144 46 : IF (trial_lambda < minimum_lambda) EXIT
145 46 : CALL set_trial_density(rho_ao, scf_env%p_mix_new, direction, trial_lambda)
146 46 : CALL evaluate_trial(qs_env, scf_env, rho, ks_env)
147 46 : evaluations = evaluations + 1
148 46 : current_lambda = trial_lambda
149 46 : trial_energy = energy%total
150 46 : CALL directional_derivative(matrix_ks, direction, gradient1)
151 46 : trial_accepted = .FALSE.
152 46 : IF (ieee_is_finite(gradient1)) THEN
153 : trial_accepted = armijo_satisfied(base_energy%total, gradient0, trial_lambda, &
154 46 : trial_energy, armijo_factor)
155 : END IF
156 46 : IF (trial_accepted) THEN
157 34 : applied = .TRUE.
158 34 : lambda = trial_lambda
159 34 : IF (gradient1 < 0.0_dp) THEN
160 26 : predictor_lambda = MIN(1.0_dp, 2.0_dp*lambda)
161 : ELSE
162 8 : predictor_lambda = MIN(1.0_dp, predictor_growth*lambda)
163 : END IF
164 34 : predictor_valid = .TRUE.
165 34 : status = status_success
166 34 : EXIT oda_search
167 : END IF
168 :
169 : ! Scale the interval to [0,1] before using the Hermite model. The model only
170 : ! proposes a smaller trial; every accepted density is still evaluated explicitly.
171 : CALL cubic_step(base_energy%total, current_lambda*gradient0, trial_energy, &
172 12 : current_lambda*gradient1, proposed_fraction, model_valid)
173 58 : IF (model_valid) THEN
174 12 : proposed_lambda = current_lambda*proposed_fraction
175 : trial_lambda = MIN(maximum_backtrack_fraction*current_lambda, &
176 12 : MAX(minimum_backtrack_fraction*current_lambda, proposed_lambda))
177 : ELSE
178 0 : trial_lambda = 0.5_dp*current_lambda
179 : END IF
180 : END DO
181 :
182 : ! Preserve the pre-ODA raw update as an explicit fallback. Restore its internally
183 : ! consistent density and Fock state after unsuccessful interior trials.
184 0 : IF (ABS(current_lambda - 1.0_dp) > 64.0_dp*EPSILON(1.0_dp)) THEN
185 0 : CALL set_trial_density(rho_ao, scf_env%p_mix_new, direction, 1.0_dp)
186 0 : CALL evaluate_trial(qs_env, scf_env, rho, ks_env)
187 0 : evaluations = evaluations + 1
188 0 : CALL directional_derivative(matrix_ks, direction, gradient1)
189 : END IF
190 0 : lambda = 1.0_dp
191 0 : status = status_no_acceptable
192 : END BLOCK oda_search
193 :
194 34 : IF (evaluations > 0) THEN
195 34 : trial_energy = energy%total
196 34 : evaluated_energy = energy
197 34 : state_evaluated = .TRUE.
198 : END IF
199 34 : IF (applied) THEN
200 : ! Keep the candidate storage aligned for consumers that inspect it before the next diagonalization.
201 610 : DO icell = 1, SIZE(rho_ao, 2)
202 1690 : DO ispin = 1, SIZE(rho_ao, 1)
203 1080 : CALL dbcsr_get_info(scf_env%p_mix_new(ispin, icell)%matrix, name=name)
204 1656 : CALL dbcsr_copy(scf_env%p_mix_new(ispin, icell)%matrix, rho_ao(ispin, icell)%matrix, name=name)
205 : END DO
206 : END DO
207 : END IF
208 34 : scf_env%oda_lambda = lambda
209 34 : scf_env%oda_energy = trial_energy
210 34 : scf_env%oda_gradient0 = gradient0
211 34 : scf_env%oda_gradient1 = gradient1
212 34 : scf_env%oda_evaluations = evaluations
213 34 : scf_env%oda_status = status
214 :
215 : ! Keep the ordinary SCF iteration's energy convention: its printed energy components
216 : ! belong to the input density P0. The accepted ODA energy is returned separately.
217 34 : energy = base_energy
218 34 : IF (ASSOCIATED(direction)) CALL dbcsr_deallocate_matrix_set(direction)
219 34 : CALL timestop(handle)
220 :
221 34 : END SUBROUTINE qs_scf_oda_apply
222 :
223 : ! **************************************************************************************************
224 : !> \brief Propose the minimum of the cubic Hermite energy model on [0,1].
225 : !> \param energy0 Energy at lambda=0.
226 : !> \param gradient0 Directional derivative at lambda=0.
227 : !> \param energy1 Energy at lambda=1.
228 : !> \param gradient1 Directional derivative at lambda=1.
229 : !> \param lambda Proposed damping parameter.
230 : !> \param valid Whether a finite model minimum was obtained.
231 : ! **************************************************************************************************
232 12 : PURE SUBROUTINE cubic_step(energy0, gradient0, energy1, gradient1, lambda, valid)
233 :
234 : REAL(KIND=dp), INTENT(IN) :: energy0, gradient0, energy1, gradient1
235 : REAL(KIND=dp), INTENT(OUT) :: lambda
236 : LOGICAL, INTENT(OUT) :: valid
237 :
238 : INTEGER :: i, nroots
239 : REAL(KIND=dp) :: a, b, best_value, discriminant, energy_delta, model_value, q, quadratic_a, &
240 : quadratic_b, quadratic_c, root, scale, sqrt_discriminant, tolerance
241 : REAL(KIND=dp), DIMENSION(2) :: roots
242 :
243 12 : lambda = 1.0_dp
244 12 : valid = .FALSE.
245 : IF (.NOT. ieee_is_finite(energy0) .OR. .NOT. ieee_is_finite(gradient0) .OR. &
246 12 : .NOT. ieee_is_finite(energy1) .OR. .NOT. ieee_is_finite(gradient1)) RETURN
247 :
248 12 : energy_delta = energy1 - energy0
249 12 : scale = MAX(1.0_dp, ABS(energy_delta), ABS(gradient0), ABS(gradient1))
250 12 : tolerance = 128.0_dp*EPSILON(1.0_dp)*scale
251 12 : IF (gradient0 >= -tolerance) RETURN
252 :
253 : ! p(lambda)-E0 = g0*lambda + a*lambda**2 + b*lambda**3
254 12 : a = 3.0_dp*energy_delta - 2.0_dp*gradient0 - gradient1
255 12 : b = gradient0 + gradient1 - 2.0_dp*energy_delta
256 12 : IF (.NOT. ieee_is_finite(a) .OR. .NOT. ieee_is_finite(b)) RETURN
257 :
258 12 : quadratic_a = 3.0_dp*b
259 12 : quadratic_b = 2.0_dp*a
260 12 : quadratic_c = gradient0
261 12 : IF (.NOT. ieee_is_finite(quadratic_a) .OR. .NOT. ieee_is_finite(quadratic_b) .OR. &
262 : .NOT. ieee_is_finite(quadratic_c)) RETURN
263 12 : roots = 0.0_dp
264 12 : nroots = 0
265 :
266 12 : IF (ABS(quadratic_a) <= tolerance) THEN
267 2 : IF (ABS(quadratic_b) > tolerance) THEN
268 2 : nroots = 1
269 2 : roots(1) = -quadratic_c/quadratic_b
270 : END IF
271 : ELSE
272 10 : discriminant = quadratic_b*quadratic_b - 4.0_dp*quadratic_a*quadratic_c
273 10 : IF (.NOT. ieee_is_finite(discriminant)) RETURN
274 10 : IF (discriminant >= -tolerance*scale) THEN
275 10 : discriminant = MAX(0.0_dp, discriminant)
276 10 : sqrt_discriminant = SQRT(discriminant)
277 10 : q = -0.5_dp*(quadratic_b + SIGN(sqrt_discriminant, quadratic_b))
278 10 : IF (ABS(q) > tolerance) THEN
279 10 : nroots = 2
280 10 : roots(1) = q/quadratic_a
281 10 : roots(2) = quadratic_c/q
282 : ELSE
283 0 : nroots = 1
284 0 : roots(1) = -quadratic_b/(2.0_dp*quadratic_a)
285 : END IF
286 : END IF
287 : END IF
288 :
289 12 : best_value = energy_delta
290 34 : DO i = 1, nroots
291 22 : root = roots(i)
292 22 : IF (.NOT. ieee_is_finite(root)) CYCLE
293 22 : IF (root <= 0.0_dp .OR. root >= 1.0_dp) CYCLE
294 12 : model_value = gradient0*root + a*root*root + b*root*root*root
295 12 : IF (.NOT. ieee_is_finite(model_value)) CYCLE
296 24 : IF (model_value < best_value - tolerance) THEN
297 12 : best_value = model_value
298 12 : lambda = root
299 : END IF
300 : END DO
301 :
302 12 : lambda = MIN(1.0_dp, MAX(0.0_dp, lambda))
303 12 : valid = .TRUE.
304 :
305 : END SUBROUTINE cubic_step
306 :
307 : ! **************************************************************************************************
308 : !> \brief Check an Armijo sufficient-decrease condition with roundoff tolerance.
309 : !> \param energy0 Energy at lambda=0.
310 : !> \param gradient0 Directional derivative at lambda=0.
311 : !> \param lambda Trial damping parameter.
312 : !> \param trial_energy Actual energy at the trial density.
313 : !> \param armijo_factor Armijo factor in (0,1).
314 : !> \return Whether the trial satisfies sufficient decrease within roundoff tolerance.
315 : ! **************************************************************************************************
316 46 : PURE LOGICAL FUNCTION armijo_satisfied(energy0, gradient0, lambda, trial_energy, &
317 : armijo_factor) RESULT(satisfied)
318 :
319 : REAL(KIND=dp), INTENT(IN) :: energy0, gradient0, lambda, &
320 : trial_energy, armijo_factor
321 :
322 : REAL(KIND=dp) :: tolerance
323 :
324 46 : satisfied = .FALSE.
325 : IF (.NOT. ieee_is_finite(energy0) .OR. .NOT. ieee_is_finite(gradient0) .OR. &
326 46 : .NOT. ieee_is_finite(lambda) .OR. .NOT. ieee_is_finite(trial_energy) .OR. &
327 : .NOT. ieee_is_finite(armijo_factor)) RETURN
328 : IF (gradient0 >= 0.0_dp .OR. lambda <= 0.0_dp .OR. lambda > 1.0_dp .OR. &
329 46 : armijo_factor <= 0.0_dp .OR. armijo_factor >= 1.0_dp) RETURN
330 :
331 46 : tolerance = 128.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(energy0), ABS(trial_energy))
332 46 : satisfied = trial_energy <= energy0 + armijo_factor*lambda*gradient0 + tolerance
333 :
334 46 : END FUNCTION armijo_satisfied
335 :
336 : ! **************************************************************************************************
337 : !> \brief Allocate and form Delta P = P1-P0.
338 : !> \param endpoint Diagonalized raw endpoint P1.
339 : !> \param base Accepted density P0.
340 : !> \param direction Allocated density direction P1-P0.
341 : ! **************************************************************************************************
342 34 : SUBROUTINE create_direction(endpoint, base, direction)
343 :
344 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: endpoint, base
345 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: direction
346 :
347 : INTEGER :: icell, ispin
348 :
349 34 : CALL dbcsr_allocate_matrix_set(direction, SIZE(endpoint, 1), SIZE(endpoint, 2))
350 610 : DO icell = 1, SIZE(endpoint, 2)
351 1690 : DO ispin = 1, SIZE(endpoint, 1)
352 1080 : ALLOCATE (direction(ispin, icell)%matrix)
353 : CALL dbcsr_create(direction(ispin, icell)%matrix, template=endpoint(ispin, icell)%matrix, &
354 1080 : name="ODA DENSITY DIRECTION")
355 1080 : CALL dbcsr_copy(direction(ispin, icell)%matrix, endpoint(ispin, icell)%matrix)
356 : CALL dbcsr_add(direction(ispin, icell)%matrix, base(ispin, icell)%matrix, &
357 1656 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
358 : END DO
359 : END DO
360 :
361 34 : END SUBROUTINE create_direction
362 :
363 : ! **************************************************************************************************
364 : !> \brief Set P(lambda) = P1 + (lambda-1) Delta P.
365 : !> \param density Trial density to overwrite.
366 : !> \param endpoint Diagonalized raw endpoint P1.
367 : !> \param direction Density direction P1-P0.
368 : !> \param lambda Trial damping parameter.
369 : ! **************************************************************************************************
370 46 : SUBROUTINE set_trial_density(density, endpoint, direction, lambda)
371 :
372 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: density
373 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: endpoint, direction
374 : REAL(KIND=dp), INTENT(IN) :: lambda
375 :
376 : CHARACTER(len=default_string_length) :: name
377 : INTEGER :: icell, ispin
378 :
379 1106 : DO icell = 1, SIZE(density, 2)
380 3154 : DO ispin = 1, SIZE(density, 1)
381 2048 : CALL dbcsr_get_info(density(ispin, icell)%matrix, name=name)
382 2048 : CALL dbcsr_copy(density(ispin, icell)%matrix, endpoint(ispin, icell)%matrix, name=name)
383 3108 : IF (lambda /= 1.0_dp) THEN
384 : CALL dbcsr_add(density(ispin, icell)%matrix, direction(ispin, icell)%matrix, &
385 2016 : alpha_scalar=1.0_dp, beta_scalar=lambda - 1.0_dp)
386 : END IF
387 : END DO
388 : END DO
389 :
390 46 : END SUBROUTINE set_trial_density
391 :
392 : ! **************************************************************************************************
393 : !> \brief Sum Tr(F Delta P) over spin channels and real-space image cells.
394 : !> \param fock Raw Fock matrices.
395 : !> \param direction Density direction P1-P0.
396 : !> \param derivative Directional energy derivative.
397 : ! **************************************************************************************************
398 80 : SUBROUTINE directional_derivative(fock, direction, derivative)
399 :
400 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: fock, direction
401 : REAL(KIND=dp), INTENT(OUT) :: derivative
402 :
403 : INTEGER :: icell, ispin
404 : REAL(KIND=dp) :: contribution
405 :
406 80 : derivative = 0.0_dp
407 : ! The k-point density transform has already folded the k-point weights into
408 : ! each real-space density cell, matching the contraction used by core energies.
409 1716 : DO icell = 1, SIZE(fock, 2)
410 4844 : DO ispin = 1, SIZE(fock, 1)
411 3128 : CALL dbcsr_dot(fock(ispin, icell)%matrix, direction(ispin, icell)%matrix, contribution)
412 4764 : derivative = derivative + contribution
413 : END DO
414 : END DO
415 :
416 80 : END SUBROUTINE directional_derivative
417 :
418 : ! **************************************************************************************************
419 : !> \brief Rebuild the density-dependent energy and raw Fock matrix for one ODA trial.
420 : !> \param qs_env Quickstep environment.
421 : !> \param scf_env SCF environment.
422 : !> \param rho Trial density.
423 : !> \param ks_env Kohn-Sham environment.
424 : ! **************************************************************************************************
425 46 : SUBROUTINE evaluate_trial(qs_env, scf_env, rho, ks_env)
426 :
427 : TYPE(qs_environment_type), POINTER :: qs_env
428 : TYPE(qs_scf_env_type), POINTER :: scf_env
429 : TYPE(qs_rho_type), POINTER :: rho
430 : TYPE(qs_ks_env_type), POINTER :: ks_env
431 :
432 46 : CALL qs_scf_rho_update(rho, qs_env, scf_env, ks_env, mix_rho=.FALSE.)
433 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE., &
434 46 : print_active=.FALSE.)
435 :
436 46 : END SUBROUTINE evaluate_trial
437 :
438 : END MODULE qs_scf_oda
|
