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 orbital transformations
10 : !> \par History
11 : !> None
12 : !> \author Joost VandeVondele (09.2002)
13 : ! **************************************************************************************************
14 : MODULE qs_ot_minimizer
15 :
16 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
17 : dbcsr_copy,&
18 : dbcsr_get_info,&
19 : dbcsr_p_type,&
20 : dbcsr_release,&
21 : dbcsr_scale,&
22 : dbcsr_set,&
23 : dbcsr_type
24 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
25 : dbcsr_dot_batch,&
26 : dbcsr_init_random
27 : USE cp_dbcsr_operations, ONLY: dbcsr_sum_many
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_unit_nr,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: high_print_level
32 : USE ieee_arithmetic, ONLY: ieee_is_finite
33 : USE kinds, ONLY: dp,&
34 : int_8
35 : USE mathlib, ONLY: diamat_all
36 : USE message_passing, ONLY: mp_comm_type,&
37 : mp_para_env_type
38 : USE preconditioner, ONLY: apply_preconditioner,&
39 : apply_preconditioner_complex_batch
40 : USE preconditioner_types, ONLY: preconditioner_p_type
41 : USE qs_ot, ONLY: qs_ot_get_derivative,&
42 : qs_ot_get_derivative_ref
43 : USE qs_ot_kp, ONLY: qs_ot_get_derivative_complex,&
44 : qs_ot_get_derivative_local_batch,&
45 : qs_ot_get_derivative_ref_complex,&
46 : qs_ot_prepare_complex_tangent_metric,&
47 : qs_ot_prepare_complex_tangent_metric_batch,&
48 : qs_ot_project_gradient
49 : USE qs_ot_types, ONLY: qs_ot_kpoint_preconditioner_scale,&
50 : qs_ot_reset_history,&
51 : qs_ot_type
52 : #include "./base/base_uses.f90"
53 :
54 : IMPLICIT NONE
55 :
56 : PRIVATE
57 :
58 : PUBLIC :: broyden_history_restart_required, &
59 : cg_history_restart_required, &
60 : lbfgs_curvature_damping_shift, &
61 : lbfgs_history_restart_required, &
62 : lbfgs_response_secant_parameters, &
63 : lbfgs_step_restart_required, &
64 : ot_mermin_response_assess, &
65 : ot_mermin_response_candidate_preferred, &
66 : ot_mermin_response_compare, &
67 : ot_mermin_response_preparation_needed, &
68 : ot_mermin_response_probe, &
69 : ot_mermin_response_shadow_followup, &
70 : ot_mermin_secant_curvature, &
71 : ot_mini, &
72 : ot_mini_prepare, &
73 : ot_mini_prepare_gradient
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_minimizer'
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief Decide whether a non-descent Broyden step invalidates its secant history.
81 : !> \param non_descent whether the proposed Broyden direction is not a descent direction
82 : !> \param forget_history explicit user request to discard an inconsistent history
83 : !> \param do_ener whether the product vector contains Mermin auxiliary-energy variables
84 : !> \return true when the Broyden history has to be restarted
85 : ! **************************************************************************************************
86 12 : PURE ELEMENTAL LOGICAL FUNCTION broyden_history_restart_required( &
87 : non_descent, forget_history, do_ener)
88 : LOGICAL, INTENT(IN) :: non_descent, forget_history, do_ener
89 :
90 12 : broyden_history_restart_required = non_descent .AND. (forget_history .OR. do_ener)
91 12 : END FUNCTION broyden_history_restart_required
92 :
93 : ! **************************************************************************************************
94 : !> \brief Decide whether an unresolved accepted energy change invalidates CG conjugacy.
95 : !>
96 : !> The occupation-preconditioned Mermin direction is a coupled Schur step rather than the action of
97 : !> one fixed positive-definite metric on the physical gradient. Once the accepted free-energy change
98 : !> is below floating-point resolution, its nonlinear conjugacy cannot be calibrated reliably. The
99 : !> current preconditioned descent direction remains valid, but the history contribution is discarded.
100 : !> \param occupation_preconditioned whether the product direction uses occupation preconditioning
101 : !> \param current_energy current Mermin free energy
102 : !> \param reference_energy Mermin free energy before the accepted line search
103 : !> \return true when the CG history contribution has to be discarded
104 : ! **************************************************************************************************
105 23136 : PURE ELEMENTAL LOGICAL FUNCTION cg_history_restart_required( &
106 : occupation_preconditioned, current_energy, reference_energy)
107 : LOGICAL, INTENT(IN) :: occupation_preconditioned
108 : REAL(KIND=dp), INTENT(IN) :: current_energy, reference_energy
109 :
110 : REAL(KIND=dp) :: energy_scale
111 :
112 23136 : energy_scale = MAX(1.0_dp, ABS(current_energy), ABS(reference_energy))
113 : cg_history_restart_required = occupation_preconditioned .AND. &
114 : ABS(current_energy - reference_energy) <= &
115 23136 : 64.0_dp*EPSILON(1.0_dp)*energy_scale
116 :
117 23136 : END FUNCTION cg_history_restart_required
118 :
119 : !
120 : ! the minimizer interface
121 : ! should present all possible modes of minimization
122 : ! these include CG SD DIIS
123 : !
124 : !
125 : ! IN the case of nspin != 1 we have a gradient that is distributed over different qs_ot_env.
126 : ! still things remain basically the same, since there are no constraints between the different qs_ot_env
127 : ! we only should take care that the various scalar products are taken over the full vectors.
128 : ! all the information needed and collected can be stored in the fist qs_ot_env only
129 : ! (indicating that the data type for the gradient/position and minization should be separated)
130 : !
131 : ! **************************************************************************************************
132 : !> \brief ...
133 : !> \param qs_ot_env ...
134 : !> \param matrix_hc ...
135 : !> \param matrix_hc_im ...
136 : !> \param matrix_hc_physical occupation-weighted derivative used for slopes and convergence
137 : !> \param matrix_hc_physical_im imaginary component of matrix_hc_physical
138 : !> \param para_env_inter_kp communicator between distributed k-point groups
139 : !> \param gradient_only return after evaluating the current OT gradient
140 : !> \param gradient_prepared reuse a gradient evaluated by ot_mini_prepare_gradient
141 : !> \param restart_history pending restart from ot_mini_prepare; absent prepares here
142 : ! **************************************************************************************************
143 97516 : SUBROUTINE ot_mini(qs_ot_env, matrix_hc, matrix_hc_im, &
144 : matrix_hc_physical, matrix_hc_physical_im, para_env_inter_kp, &
145 : gradient_only, gradient_prepared, restart_history)
146 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
147 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
148 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
149 : POINTER :: matrix_hc_im, matrix_hc_physical, &
150 : matrix_hc_physical_im
151 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
152 : LOGICAL, INTENT(IN), OPTIONAL :: gradient_only, gradient_prepared, &
153 : restart_history
154 :
155 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_mini'
156 :
157 : INTEGER :: handle, ispin, npasses, nspin, pass
158 : LOGICAL :: do_ener, do_ks, evaluate_gradient_only, &
159 : local_gradient, restart, reuse_gradient
160 : REAL(KIND=dp) :: tmp
161 : TYPE(dbcsr_type), POINTER :: hc_im, hc_re
162 :
163 97516 : CALL timeset(routineN, handle)
164 :
165 97516 : evaluate_gradient_only = .FALSE.
166 97516 : reuse_gradient = .FALSE.
167 97516 : IF (PRESENT(gradient_only)) evaluate_gradient_only = gradient_only
168 97516 : IF (PRESENT(gradient_prepared)) reuse_gradient = gradient_prepared
169 3050 : CPASSERT(.NOT. (evaluate_gradient_only .AND. reuse_gradient))
170 :
171 97516 : nspin = SIZE(qs_ot_env)
172 :
173 97516 : do_ks = qs_ot_env(1)%settings%ks
174 97516 : do_ener = qs_ot_env(1)%settings%do_ener
175 :
176 97516 : IF (PRESENT(restart_history)) THEN
177 84374 : restart = restart_history
178 : ELSE
179 13142 : CALL ot_mini_prepare(qs_ot_env, restart, para_env_inter_kp)
180 : END IF
181 97516 : IF (restart) THEN
182 17234 : DO ispin = 1, nspin
183 17234 : CALL qs_ot_reset_history(qs_ot_env(ispin))
184 : END DO
185 : reuse_gradient = .FALSE.
186 : END IF
187 :
188 97516 : npasses = 1
189 97516 : IF (PRESENT(matrix_hc_physical)) THEN
190 992 : CPASSERT(qs_ot_env(1)%settings%occupation_preconditioner)
191 992 : CPASSERT(SIZE(matrix_hc_physical) == nspin)
192 : npasses = 2
193 : END IF
194 :
195 97516 : qs_ot_env(1)%OT_METHOD_FULL = ""
196 :
197 : ! compute the gradient for the variables x
198 97516 : IF (.NOT. reuse_gradient .AND. .NOT. qs_ot_env(1)%energy_only) THEN
199 74970 : qs_ot_env(1)%gradient = 0.0_dp
200 74970 : local_gradient = .FALSE.
201 74970 : IF (do_ks .AND. qs_ot_env(1)%has_complex_kpoint_state) THEN
202 2410 : IF (.NOT. PRESENT(matrix_hc_im)) THEN
203 0 : CPABORT("Complex k-point OT derivative requires imaginary H(k)*C(k).")
204 : END IF
205 2410 : local_gradient = qs_ot_env(1)%para_env%num_pe == 1
206 2410 : IF (local_gradient) THEN
207 : CALL qs_ot_get_derivative_local_batch(qs_ot_env, matrix_hc, matrix_hc_im, &
208 1256 : matrix_hc_physical, matrix_hc_physical_im)
209 : END IF
210 : END IF
211 162702 : DO ispin = 1, nspin
212 87732 : IF (do_ks .AND. .NOT. local_gradient) THEN
213 170732 : DO pass = 1, npasses
214 85671 : hc_re => matrix_hc(ispin)%matrix
215 85671 : IF (pass == 2) hc_re => matrix_hc_physical(ispin)%matrix
216 85671 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
217 2386 : CPASSERT(PRESENT(matrix_hc_im))
218 2386 : hc_im => matrix_hc_im(ispin)%matrix
219 2386 : IF (pass == 2) THEN
220 610 : CPASSERT(PRESENT(matrix_hc_physical_im))
221 610 : hc_im => matrix_hc_physical_im(ispin)%matrix
222 : END IF
223 3690 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
224 : CASE ("TOD")
225 1304 : CALL qs_ot_get_derivative_complex(hc_re, hc_im, qs_ot_env(ispin))
226 : CASE ("REF")
227 1082 : CALL qs_ot_get_derivative_ref_complex(hc_re, hc_im, qs_ot_env(ispin))
228 : CASE DEFAULT
229 2386 : CPABORT("Complex k-point OT derivative requires ALGORITHM STRICT or IRAC")
230 : END SELECT
231 : ELSE
232 163038 : SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
233 : CASE ("TOD")
234 : CALL qs_ot_get_derivative(hc_re, qs_ot_env(ispin)%matrix_x, &
235 : qs_ot_env(ispin)%matrix_sx, &
236 79753 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
237 : CASE ("REF")
238 : CALL qs_ot_get_derivative_ref(hc_re, &
239 : qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, &
240 3532 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
241 : CASE DEFAULT
242 83285 : CPABORT("ALGORITHM NYI")
243 : END SELECT
244 : END IF
245 170732 : IF (pass < npasses) THEN
246 : ! Save the search gradient; the final pass supplies physical X/rotation gradients.
247 610 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%matrix_preconditioned_gx))
248 610 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_preconditioned_gx, qs_ot_env(ispin)%matrix_gx)
249 610 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
250 610 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%matrix_preconditioned_gx_im))
251 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_preconditioned_gx_im, &
252 610 : qs_ot_env(ispin)%matrix_gx_im)
253 : END IF
254 : END IF
255 : END DO
256 : END IF
257 : ! and also the gradient along the direction
258 162702 : IF (qs_ot_env(1)%use_dx) THEN
259 82764 : IF (do_ks) THEN
260 82764 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
261 82764 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
262 82764 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
263 4395 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_dx_im, tmp)
264 4395 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
265 : END IF
266 82764 : IF (qs_ot_env(1)%settings%do_rotation) THEN
267 2969 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
268 2969 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + 0.5_dp*tmp
269 2969 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
270 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
271 1373 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
272 1373 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + 0.5_dp*tmp
273 : END IF
274 : END IF
275 : END IF
276 82764 : IF (do_ener) THEN
277 8322 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
278 966 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
279 : END IF
280 : ELSE
281 4968 : IF (do_ks) THEN
282 4968 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
283 4968 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
284 4968 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
285 52 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_gx_im, tmp)
286 52 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
287 : END IF
288 4968 : IF (qs_ot_env(1)%settings%do_rotation) THEN
289 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
290 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - 0.5_dp*tmp
291 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
292 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
293 0 : qs_ot_env(ispin)%rot_mat_gx_im, tmp)
294 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - 0.5_dp*tmp
295 : END IF
296 : END IF
297 : END IF
298 4968 : IF (do_ener) THEN
299 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
300 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
301 : END IF
302 : END IF
303 : END DO
304 74970 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%gradient)
305 : END IF
306 :
307 97516 : IF (evaluate_gradient_only) THEN
308 648 : CALL timestop(handle)
309 648 : RETURN
310 : END IF
311 :
312 141306 : SELECT CASE (qs_ot_env(1)%settings%OT_METHOD)
313 : CASE ("CG")
314 44438 : IF (current_point_is_fine(qs_ot_env)) THEN
315 23132 : qs_ot_env(1)%OT_METHOD_FULL = "OT CG"
316 23132 : CALL ot_new_cg_direction(qs_ot_env, para_env_inter_kp)
317 23132 : qs_ot_env(1)%line_search_count = 0
318 : ELSE
319 21306 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
320 : END IF
321 44438 : CALL do_line_search(qs_ot_env)
322 : CASE ("SD")
323 164 : IF (current_point_is_fine(qs_ot_env)) THEN
324 50 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
325 50 : CALL ot_new_sd_direction(qs_ot_env, para_env_inter_kp)
326 50 : qs_ot_env(1)%line_search_count = 0
327 : ELSE
328 114 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
329 : END IF
330 164 : CALL do_line_search(qs_ot_env)
331 : CASE ("DIIS")
332 51136 : qs_ot_env(1)%OT_METHOD_FULL = "OT DIIS"
333 51136 : CALL ot_diis_step(qs_ot_env, para_env_inter_kp)
334 : CASE ("BROY")
335 192 : qs_ot_env(1)%OT_METHOD_FULL = "OT BROY"
336 192 : CALL ot_broyden_step(qs_ot_env, para_env_inter_kp)
337 : CASE ("LBFG")
338 938 : IF (current_point_is_fine(qs_ot_env)) THEN
339 460 : qs_ot_env(1)%OT_METHOD_FULL = "OT LBFGS"
340 460 : CALL ot_new_lbfgs_direction(qs_ot_env, para_env_inter_kp)
341 460 : qs_ot_env(1)%line_search_count = 0
342 : ELSE
343 478 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
344 : END IF
345 938 : CALL do_line_search(qs_ot_env)
346 : CASE DEFAULT
347 96868 : CPABORT("OT_METHOD NYI")
348 : END SELECT
349 :
350 96868 : CALL timestop(handle)
351 :
352 : END SUBROUTINE ot_mini
353 :
354 : ! **************************************************************************************************
355 : !> \brief Prepare projection metrics and synchronize pending history restarts before building gradients.
356 : !> \param qs_ot_env OT channel environments
357 : !> \param restart_history whether ot_mini must restart and evaluate a gradient, even on an energy-only step
358 : !> \param para_env_inter_kp communicator between distributed k-point groups
359 : ! **************************************************************************************************
360 96868 : SUBROUTINE ot_mini_prepare(qs_ot_env, restart_history, para_env_inter_kp)
361 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
362 : LOGICAL, INTENT(OUT) :: restart_history
363 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
364 :
365 : INTEGER :: ispin, restart
366 :
367 96868 : restart = 0
368 96868 : CALL qs_ot_prepare_complex_tangent_metric_batch(qs_ot_env)
369 211388 : DO ispin = 1, SIZE(qs_ot_env)
370 114520 : IF (qs_ot_env(1)%settings%ks .AND. qs_ot_env(ispin)%settings%ot_algorithm == 'TOD' .AND. &
371 : qs_ot_env(ispin)%has_complex_kpoint_state) THEN
372 3314 : CALL qs_ot_prepare_complex_tangent_metric(qs_ot_env(ispin))
373 : END IF
374 211388 : IF (qs_ot_env(ispin)%preconditioner_changed) restart = 1
375 : END DO
376 96868 : IF (PRESENT(para_env_inter_kp)) THEN
377 3050 : IF (ASSOCIATED(para_env_inter_kp)) CALL para_env_inter_kp%max(restart)
378 : END IF
379 96868 : restart_history = restart /= 0
380 96868 : END SUBROUTINE ot_mini_prepare
381 :
382 : ! **************************************************************************************************
383 : !> \brief Evaluate the current OT derivative without advancing the minimizer.
384 : !> \param qs_ot_env OT channel environments
385 : !> \param matrix_hc real H*C products
386 : !> \param matrix_hc_im imaginary H*C products
387 : !> \param matrix_hc_physical occupation-weighted real H*C products
388 : !> \param matrix_hc_physical_im occupation-weighted imaginary H*C products
389 : !> \param para_env_inter_kp communicator between distributed K-point groups
390 : !> \param restart_history pending restart from ot_mini_prepare; absent prepares here
391 : ! **************************************************************************************************
392 648 : SUBROUTINE ot_mini_prepare_gradient(qs_ot_env, matrix_hc, matrix_hc_im, &
393 : matrix_hc_physical, matrix_hc_physical_im, para_env_inter_kp, &
394 : restart_history)
395 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
396 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
397 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
398 : POINTER :: matrix_hc_im, matrix_hc_physical, &
399 : matrix_hc_physical_im
400 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
401 : LOGICAL, INTENT(IN), OPTIONAL :: restart_history
402 :
403 : CALL ot_mini(qs_ot_env, matrix_hc, matrix_hc_im=matrix_hc_im, &
404 : matrix_hc_physical=matrix_hc_physical, &
405 : matrix_hc_physical_im=matrix_hc_physical_im, &
406 648 : para_env_inter_kp=para_env_inter_kp, gradient_only=.TRUE., restart_history=restart_history)
407 :
408 648 : END SUBROUTINE ot_mini_prepare_gradient
409 :
410 : ! **************************************************************************************************
411 : !> \brief Apply a channel's orbital preconditioner, or a scaled identity when absent.
412 : !> \param qs_ot_env current OT channel
413 : !> \param matrix_in input real part
414 : !> \param matrix_out output real part
415 : !> \param matrix_in_im input imaginary part, unused for real channels
416 : !> \param matrix_out_im output imaginary part, unused for real channels
417 : !> \param identity_scale scale used only when the preconditioner is absent
418 : ! **************************************************************************************************
419 87294 : SUBROUTINE ot_apply_channel_preconditioner(qs_ot_env, matrix_in, matrix_out, &
420 : matrix_in_im, matrix_out_im, identity_scale)
421 : TYPE(qs_ot_type), INTENT(IN) :: qs_ot_env
422 : TYPE(dbcsr_type), POINTER :: matrix_in, matrix_out, matrix_in_im, &
423 : matrix_out_im
424 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: identity_scale
425 :
426 : REAL(KIND=dp) :: scale
427 :
428 87294 : scale = 1.0_dp
429 87294 : IF (ASSOCIATED(qs_ot_env%preconditioner)) THEN
430 76696 : IF (qs_ot_env%has_complex_kpoint_state) THEN
431 3909 : CALL apply_preconditioner(qs_ot_env%preconditioner, matrix_in, matrix_in_im, matrix_out, matrix_out_im)
432 3909 : CPASSERT(qs_ot_env%kpoint_weight > 0.0_dp)
433 3909 : scale = qs_ot_kpoint_preconditioner_scale(qs_ot_env%kpoint_weight)
434 : ELSE
435 72787 : CALL apply_preconditioner(qs_ot_env%preconditioner, matrix_in, matrix_out)
436 : END IF
437 : ELSE
438 10598 : CALL dbcsr_copy(matrix_out, matrix_in)
439 10598 : IF (qs_ot_env%has_complex_kpoint_state) CALL dbcsr_copy(matrix_out_im, matrix_in_im)
440 10598 : IF (PRESENT(identity_scale)) scale = identity_scale
441 : END IF
442 87294 : IF (scale /= 1.0_dp) THEN
443 2205 : CALL dbcsr_scale(matrix_out, scale)
444 2205 : IF (qs_ot_env%has_complex_kpoint_state) CALL dbcsr_scale(matrix_out_im, scale)
445 : END IF
446 87294 : END SUBROUTINE ot_apply_channel_preconditioner
447 :
448 : ! **************************************************************************************************
449 : !> \brief Precondition orbital gradients across local channels, retaining the distributed fallback.
450 : !> \param qs_ot_env OT channels
451 : !> \param matrix_out output pointers, indexed by real/imaginary component and channel
452 : !> \param use_occupation_response override the channel's occupation-response setting
453 : ! **************************************************************************************************
454 74502 : SUBROUTINE ot_precondition_gradients(qs_ot_env, matrix_out, use_occupation_response)
455 : TYPE(qs_ot_type), INTENT(IN) :: qs_ot_env(:)
456 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_out(:, :)
457 : LOGICAL, INTENT(IN), OPTIONAL :: use_occupation_response
458 :
459 : INTEGER :: i
460 : LOGICAL :: handled, local, use_response
461 161536 : TYPE(preconditioner_p_type) :: preconditioners(SIZE(qs_ot_env))
462 335604 : TYPE(dbcsr_p_type) :: matrix_in(2, SIZE(qs_ot_env))
463 149004 : REAL(dp) :: scales(SIZE(qs_ot_env))
464 :
465 161536 : local = .TRUE.
466 161536 : scales(:) = 1.0_dp
467 161536 : DO i = 1, SIZE(qs_ot_env)
468 87034 : use_response = qs_ot_env(i)%settings%occupation_preconditioner
469 87034 : IF (PRESENT(use_occupation_response)) use_response = use_occupation_response
470 87034 : IF (use_response) THEN
471 522 : matrix_in(1, i)%matrix => qs_ot_env(i)%matrix_preconditioned_gx
472 522 : matrix_in(2, i)%matrix => qs_ot_env(i)%matrix_preconditioned_gx_im
473 : ELSE
474 86512 : matrix_in(1, i)%matrix => qs_ot_env(i)%matrix_gx
475 86512 : matrix_in(2, i)%matrix => qs_ot_env(i)%matrix_gx_im
476 : END IF
477 87034 : preconditioners(i)%preconditioner => qs_ot_env(i)%preconditioner
478 161536 : IF (qs_ot_env(i)%has_complex_kpoint_state) THEN
479 3761 : IF (ASSOCIATED(qs_ot_env(i)%preconditioner)) THEN
480 3705 : CPASSERT(qs_ot_env(i)%kpoint_weight > 0.0_dp)
481 3705 : scales(i) = qs_ot_kpoint_preconditioner_scale(qs_ot_env(i)%kpoint_weight)
482 : END IF
483 : ELSE
484 : local = .FALSE.
485 : END IF
486 : END DO
487 74502 : handled = .FALSE.
488 74502 : IF (local) CALL apply_preconditioner_complex_batch(preconditioners, matrix_in, matrix_out, handled, scales)
489 74502 : IF (handled) RETURN
490 157969 : DO i = 1, SIZE(qs_ot_env)
491 : CALL ot_apply_channel_preconditioner(qs_ot_env(i), matrix_in(1, i)%matrix, &
492 : matrix_out(1, i)%matrix, matrix_in(2, i)%matrix, &
493 157969 : matrix_out(2, i)%matrix)
494 : END DO
495 : END SUBROUTINE ot_precondition_gradients
496 :
497 : ! **************************************************************************************************
498 : !> \brief Sum a minimizer scalar over distributed k-point groups when requested.
499 : !> \param para_env_inter_kp communicator between distributed k-point groups
500 : !> \param value scalar to sum
501 : ! **************************************************************************************************
502 650111 : SUBROUTINE ot_mini_sum_kpoint_scalar(para_env_inter_kp, value)
503 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
504 : REAL(KIND=dp), INTENT(INOUT) :: value
505 :
506 650111 : IF (PRESENT(para_env_inter_kp)) THEN
507 21490 : IF (ASSOCIATED(para_env_inter_kp)) CALL para_env_inter_kp%sum(value)
508 : END IF
509 :
510 650111 : END SUBROUTINE ot_mini_sum_kpoint_scalar
511 :
512 : ! **************************************************************************************************
513 : !> \brief Normalize the convergence product on an orbital/energy scale.
514 : !> \param qs_ot_env local spin/k-point channels; updates delta in the first channel
515 : !> \param nvariables global coordinate count, weighted for complex k points
516 : !> \param para_env_inter_kp communicator between distributed k-point groups
517 : !> \param residual_product optional current, globally reduced, unscaled gradient/preconditioned-gradient product
518 : !> \note Call before overwriting the physical gradients for complex k points.
519 : ! **************************************************************************************************
520 74970 : SUBROUTINE ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp, residual_product)
521 : TYPE(qs_ot_type), DIMENSION(:), INTENT(INOUT) :: qs_ot_env
522 : REAL(KIND=dp), INTENT(OUT) :: nvariables
523 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
524 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: residual_product
525 :
526 : INTEGER :: ispin, k, n
527 : INTEGER(KIND=int_8) :: nlocal, norbital
528 : REAL(KIND=dp) :: gnorm, weight
529 :
530 74970 : gnorm = qs_ot_env(1)%gnorm
531 74970 : IF (qs_ot_env(1)%has_complex_kpoint_state) THEN
532 2410 : gnorm = ot_kpoint_convergence(qs_ot_env, para_env_inter_kp, residual_product)
533 : END IF
534 :
535 : ! Rotations contribute to the product, but do not add coordinates to the
536 : ! orbital/energy normalization and dilute the stopping criterion.
537 74970 : nvariables = 0.0_dp
538 162702 : DO ispin = 1, SIZE(qs_ot_env)
539 87732 : nlocal = 0_int_8
540 87732 : IF (qs_ot_env(1)%settings%ks) THEN
541 87732 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullrows_total=n, nfullcols_total=k)
542 87732 : norbital = INT(n, KIND=int_8)*INT(k, KIND=int_8)
543 87732 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) norbital = 2_int_8*norbital
544 : nlocal = nlocal + norbital
545 : END IF
546 87732 : IF (qs_ot_env(1)%settings%do_ener) THEN
547 966 : nlocal = nlocal + SIZE(qs_ot_env(ispin)%ener_x, KIND=int_8)
548 : END IF
549 87732 : weight = 1.0_dp
550 87732 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) weight = qs_ot_env(ispin)%kpoint_weight
551 162702 : nvariables = nvariables + weight*REAL(nlocal, KIND=dp)
552 : END DO
553 74970 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, nvariables)
554 :
555 74970 : qs_ot_env(1)%delta = 0.0_dp
556 74970 : IF (nvariables > 0.0_dp) qs_ot_env(1)%delta = SQRT(ABS(gnorm)/nvariables)
557 74970 : END SUBROUTINE ot_update_delta
558 :
559 : ! **************************************************************************************************
560 : !> \brief Evaluate the coupled residual with one k-point weight per contribution.
561 : !> \param qs_ot_env complex k-point channels with their current physical gradients
562 : !> \param para_env_inter_kp communicator between distributed k-point groups
563 : !> \param residual_product optional current, globally reduced, unscaled residual product
564 : !> \return weighted convergence product, independent of search directions and history
565 : ! **************************************************************************************************
566 2410 : FUNCTION ot_kpoint_convergence(qs_ot_env, para_env_inter_kp, residual_product) RESULT(gnorm)
567 : TYPE(qs_ot_type), DIMENSION(:), INTENT(INOUT) :: qs_ot_env
568 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
569 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: residual_product
570 : REAL(KIND=dp) :: gnorm
571 :
572 : INTEGER, PARAMETER :: identity_residual = 3, &
573 : occupation_response = 1, &
574 : physical_residual = 2
575 :
576 : INTEGER :: ispin, residual_kind, reuse_product
577 : REAL(KIND=dp) :: inverse_weight, scale, tmp
578 : TYPE(dbcsr_type), POINTER :: input_im, input_re, residual_im, &
579 : residual_re
580 : TYPE(dbcsr_type), TARGET :: identity_im, identity_re
581 :
582 2410 : IF (PRESENT(residual_product)) THEN
583 1950 : reuse_product = 0
584 1950 : IF (ieee_is_finite(residual_product)) THEN
585 1950 : IF (residual_product > 0.0_dp) reuse_product = 1
586 : END IF
587 : ! Without an attached preconditioner, the history product lacks the inverse k-point weight.
588 5703 : DO ispin = 1, SIZE(qs_ot_env)
589 : IF (.NOT. ASSOCIATED(qs_ot_env(ispin)%preconditioner) .OR. &
590 : qs_ot_env(ispin)%settings%do_ener .OR. &
591 3753 : qs_ot_env(ispin)%settings%do_rotation .OR. &
592 2813 : qs_ot_env(ispin)%settings%occupation_preconditioner) reuse_product = 0
593 : END DO
594 : ! All groups must take the same branch before the fallback's collective reductions.
595 1950 : IF (PRESENT(para_env_inter_kp)) THEN
596 1462 : IF (ASSOCIATED(para_env_inter_kp)) CALL para_env_inter_kp%min(reuse_product)
597 : END IF
598 1950 : IF (reuse_product == 1) THEN
599 : gnorm = residual_product
600 : RETURN
601 : END IF
602 : END IF
603 :
604 984 : DO residual_kind = MERGE(occupation_response, physical_residual, &
605 978 : qs_ot_env(1)%settings%occupation_preconditioner), identity_residual
606 984 : gnorm = 0.0_dp
607 2551 : DO ispin = 1, SIZE(qs_ot_env)
608 1567 : CPASSERT(qs_ot_env(ispin)%kpoint_weight > 0.0_dp)
609 1567 : inverse_weight = 1.0_dp/qs_ot_env(ispin)%kpoint_weight
610 1567 : IF (qs_ot_env(1)%settings%ks) THEN
611 1567 : input_re => qs_ot_env(ispin)%matrix_gx
612 1567 : input_im => qs_ot_env(ispin)%matrix_gx_im
613 1567 : IF (residual_kind == occupation_response) THEN
614 726 : input_re => qs_ot_env(ispin)%matrix_preconditioned_gx
615 726 : input_im => qs_ot_env(ispin)%matrix_preconditioned_gx_im
616 : END IF
617 1567 : scale = inverse_weight
618 1567 : IF (residual_kind /= identity_residual) THEN
619 : ! Derivative workspaces are free here; preserve minimizer history and directions.
620 1567 : IF (qs_ot_env(ispin)%settings%ot_algorithm == 'TOD') THEN
621 572 : residual_re => qs_ot_env(ispin)%matrix_buf_nk
622 572 : residual_im => qs_ot_env(ispin)%matrix_buf_nk_im
623 : ELSE
624 995 : residual_re => qs_ot_env(ispin)%buf1_n_k
625 995 : residual_im => qs_ot_env(ispin)%buf1_n_k_dp
626 : END IF
627 : CALL ot_apply_channel_preconditioner(qs_ot_env(ispin), input_re, residual_re, &
628 1567 : input_im, residual_im)
629 : ! The attached inverse Hessian already includes 1/w; an identity map does not.
630 1567 : IF (ASSOCIATED(qs_ot_env(ispin)%preconditioner)) scale = 1.0_dp
631 : ELSE
632 : ! STRICT projection uses the derivative buffers, so this rare fallback needs
633 : ! separate targets. Only one channel's temporary matrices are live at a time.
634 0 : CALL dbcsr_copy(identity_re, input_re)
635 0 : CALL dbcsr_copy(identity_im, input_im)
636 0 : residual_re => identity_re
637 0 : residual_im => identity_im
638 0 : IF (qs_ot_env(ispin)%settings%ot_algorithm == 'TOD') THEN
639 : CALL qs_ot_project_gradient(qs_ot_env(ispin), residual_re, residual_im, &
640 0 : use_preconditioner=.FALSE.)
641 : END IF
642 : END IF
643 1567 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, residual_re, tmp)
644 1567 : gnorm = gnorm + scale*tmp
645 1567 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, residual_im, tmp)
646 1567 : gnorm = gnorm + scale*tmp
647 1567 : IF (residual_kind == identity_residual) THEN
648 0 : CALL dbcsr_release(identity_re)
649 0 : CALL dbcsr_release(identity_im)
650 : END IF
651 1567 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
652 1383 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
653 1383 : gnorm = gnorm + 0.5_dp*inverse_weight*tmp
654 1383 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, qs_ot_env(ispin)%rot_mat_gx_im, tmp)
655 1383 : gnorm = gnorm + 0.5_dp*inverse_weight*tmp
656 : END IF
657 : END IF
658 2551 : IF (qs_ot_env(1)%settings%do_ener) THEN
659 976 : IF (residual_kind == occupation_response) THEN
660 5442 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_preconditioned_gx)
661 : ELSE
662 2922 : tmp = inverse_weight*DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
663 : END IF
664 976 : gnorm = gnorm + tmp
665 : END IF
666 : END DO
667 984 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, gnorm)
668 984 : IF (ieee_is_finite(gnorm) .AND. gnorm > 0.0_dp) EXIT
669 : ! Decide the diagnostic fallback from the complete weighted product. The minimizers
670 : ! independently test their actual directions with the original, unscaled inner products.
671 : END DO
672 3388 : END FUNCTION ot_kpoint_convergence
673 :
674 : ! **************************************************************************************************
675 : !> \brief Assess a finite-response candidate at its accepted Mermin endpoint.
676 : !> \param reference_energy Mermin energy at the start of the search
677 : !> \param current_energy Mermin energy at the accepted endpoint
678 : !> \param reference_residual physical residual at the start of the search
679 : !> \param current_residual physical residual at the accepted endpoint
680 : !> \param predicted_slope positive predicted decrease per unit line-search position
681 : !> \param predicted_curvature signed finite-response curvature
682 : !> \param position accepted line-search position
683 : !> \param default_step configured initial line-search position
684 : !> \param predicted_drop finite quadratic-model decrease
685 : !> \param measured_drop measured Mermin decrease
686 : !> \param quality measured over predicted decrease
687 : !> \param residual_ratio endpoint over reference residual
688 : !> \param good whether this is a resolved, useful accepted response sample
689 : ! **************************************************************************************************
690 20 : PURE SUBROUTINE ot_mermin_response_assess( &
691 : reference_energy, current_energy, reference_residual, current_residual, &
692 : predicted_slope, predicted_curvature, position, default_step, &
693 : predicted_drop, measured_drop, quality, residual_ratio, good)
694 : REAL(KIND=dp), INTENT(IN) :: reference_energy, current_energy, reference_residual, &
695 : current_residual, predicted_slope, predicted_curvature, position, default_step
696 : REAL(KIND=dp), INTENT(OUT) :: predicted_drop, measured_drop, quality, &
697 : residual_ratio
698 : LOGICAL, INTENT(OUT) :: good
699 :
700 : REAL(KIND=dp) :: accepted_position, resolution, scale
701 :
702 20 : predicted_drop = 0.0_dp
703 20 : measured_drop = 0.0_dp
704 20 : quality = 0.0_dp
705 20 : residual_ratio = HUGE(1.0_dp)
706 20 : good = .FALSE.
707 : IF (.NOT. ieee_is_finite(reference_energy) .OR. .NOT. ieee_is_finite(current_energy) .OR. &
708 : .NOT. ieee_is_finite(reference_residual) .OR. .NOT. ieee_is_finite(current_residual) .OR. &
709 : .NOT. ieee_is_finite(predicted_slope) .OR. .NOT. ieee_is_finite(predicted_curvature) .OR. &
710 20 : .NOT. ieee_is_finite(position) .OR. .NOT. ieee_is_finite(default_step)) RETURN
711 :
712 20 : accepted_position = ABS(position)
713 20 : measured_drop = reference_energy - current_energy
714 : predicted_drop = accepted_position*predicted_slope - &
715 20 : 0.5_dp*accepted_position**2*predicted_curvature
716 : scale = MAX(1.0_dp, ABS(reference_energy), ABS(current_energy), &
717 20 : ABS(predicted_drop), ABS(measured_drop))
718 20 : resolution = 256.0_dp*EPSILON(1.0_dp)*scale
719 20 : IF (reference_residual > TINY(reference_residual)) THEN
720 20 : residual_ratio = MAX(0.0_dp, current_residual)/reference_residual
721 : END IF
722 20 : IF (predicted_drop > resolution) quality = measured_drop/predicted_drop
723 :
724 : good = predicted_drop > resolution .AND. measured_drop > resolution .AND. &
725 : accepted_position >= 1.0E-6_dp*MAX(ABS(default_step), TINY(default_step)) .AND. &
726 26 : quality >= 0.05_dp .AND. quality <= 5.0_dp .AND. residual_ratio <= 1.0_dp
727 :
728 : END SUBROUTINE ot_mermin_response_assess
729 :
730 : ! **************************************************************************************************
731 : !> \brief Recover the total finite Mermin curvature of an accepted line-search secant.
732 : !>
733 : !> For F(alpha)=F(0)-g*alpha+0.5*kappa*alpha**2, the accepted energy difference fixes
734 : !> kappa without separating REF, rotation, occupation, and self-consistent Hxc terms.
735 : !> \param reference_energy Mermin energy at the start of the accepted search
736 : !> \param current_energy Mermin energy at its accepted endpoint
737 : !> \param predicted_slope positive directional decrease at the search origin
738 : !> \param position accepted line-search position
739 : !> \param curvature recovered total finite curvature
740 : !> \param valid whether a resolved finite curvature was recovered
741 : ! **************************************************************************************************
742 138 : PURE SUBROUTINE ot_mermin_secant_curvature( &
743 : reference_energy, current_energy, predicted_slope, position, curvature, valid)
744 : REAL(KIND=dp), INTENT(IN) :: reference_energy, current_energy, &
745 : predicted_slope, position
746 : REAL(KIND=dp), INTENT(OUT) :: curvature
747 : LOGICAL, INTENT(OUT) :: valid
748 :
749 : REAL(KIND=dp) :: accepted_position, measured_drop
750 :
751 138 : curvature = 0.0_dp
752 138 : valid = .FALSE.
753 : IF (.NOT. ieee_is_finite(reference_energy) .OR. .NOT. ieee_is_finite(current_energy) .OR. &
754 138 : .NOT. ieee_is_finite(predicted_slope) .OR. .NOT. ieee_is_finite(position)) RETURN
755 138 : accepted_position = ABS(position)
756 138 : IF (accepted_position <= SQRT(EPSILON(1.0_dp))) RETURN
757 :
758 138 : measured_drop = reference_energy - current_energy
759 138 : curvature = 2.0_dp*(accepted_position*predicted_slope - measured_drop)/accepted_position**2
760 138 : valid = ieee_is_finite(curvature)
761 138 : IF (.NOT. valid) curvature = 0.0_dp
762 :
763 : END SUBROUTINE ot_mermin_secant_curvature
764 :
765 : ! **************************************************************************************************
766 : !> \brief Compare an endpoint response shadow with the linear accepted-step model.
767 : !> \param reference_energy Mermin energy at the start of the search
768 : !> \param current_energy Mermin energy at the accepted endpoint
769 : !> \param reference_residual physical residual at the start of the search
770 : !> \param current_residual physical residual at the accepted endpoint
771 : !> \param predicted_slope positive reference directional decrease
772 : !> \param shadow_curvature endpoint finite-response curvature along the accepted direction
773 : !> \param position accepted line-search position
774 : !> \param advantage reduction of symmetric prediction error relative to the linear model
775 : !> \param residual_ratio endpoint over reference residual
776 : !> \param good whether the passive response prediction is both better and residual-consistent
777 : ! **************************************************************************************************
778 144 : PURE SUBROUTINE ot_mermin_response_compare( &
779 : reference_energy, current_energy, reference_residual, current_residual, &
780 : predicted_slope, shadow_curvature, position, advantage, residual_ratio, good)
781 : REAL(KIND=dp), INTENT(IN) :: reference_energy, current_energy, reference_residual, &
782 : current_residual, predicted_slope, shadow_curvature, position
783 : REAL(KIND=dp), INTENT(OUT) :: advantage, residual_ratio
784 : LOGICAL, INTENT(OUT) :: good
785 :
786 : REAL(KIND=dp) :: accepted_position, baseline_drop, baseline_error, candidate_drop, &
787 : candidate_error, measured_drop, resolution, scale
788 :
789 144 : advantage = 0.0_dp
790 144 : residual_ratio = HUGE(1.0_dp)
791 144 : good = .FALSE.
792 : IF (.NOT. ieee_is_finite(reference_energy) .OR. .NOT. ieee_is_finite(current_energy) .OR. &
793 : .NOT. ieee_is_finite(reference_residual) .OR. .NOT. ieee_is_finite(current_residual) .OR. &
794 144 : .NOT. ieee_is_finite(predicted_slope) .OR. .NOT. ieee_is_finite(shadow_curvature) .OR. &
795 : .NOT. ieee_is_finite(position)) RETURN
796 :
797 144 : accepted_position = ABS(position)
798 144 : measured_drop = reference_energy - current_energy
799 144 : baseline_drop = accepted_position*predicted_slope
800 144 : candidate_drop = baseline_drop - 0.5_dp*accepted_position**2*shadow_curvature
801 : scale = MAX(1.0_dp, ABS(reference_energy), ABS(current_energy), ABS(measured_drop), &
802 144 : ABS(baseline_drop), ABS(candidate_drop))
803 144 : resolution = 256.0_dp*EPSILON(1.0_dp)*scale
804 144 : IF (reference_residual > TINY(reference_residual)) THEN
805 144 : residual_ratio = MAX(0.0_dp, current_residual)/reference_residual
806 : END IF
807 144 : IF (measured_drop <= resolution .OR. baseline_drop <= resolution .OR. &
808 : candidate_drop <= resolution) RETURN
809 :
810 144 : baseline_error = ABS(measured_drop - baseline_drop)/scale
811 144 : candidate_error = ABS(measured_drop - candidate_drop)/scale
812 144 : advantage = baseline_error - candidate_error
813 160 : good = advantage > 64.0_dp*EPSILON(1.0_dp)*scale .AND. residual_ratio <= 1.05_dp
814 :
815 : END SUBROUTINE ot_mermin_response_compare
816 :
817 : ! **************************************************************************************************
818 : !> \brief Compare response and conventional directions in one accepted Mermin model.
819 : !>
820 : !> The most recent accepted conventional secant supplies the dimensionless local curvature
821 : !> ratio kappa/g. Applying that ratio to the current conventional slope transfers the
822 : !> measured line-search model without assuming that the two directions have equal norms.
823 : !> Both candidates are then evaluated at the same line-search position.
824 : !> \param baseline_slope positive current conventional directional decrease
825 : !> \param response_slope positive current response directional decrease
826 : !> \param response_curvature current finite-response curvature
827 : !> \param accepted_slope positive slope of the preceding accepted conventional search
828 : !> \param accepted_curvature measured total curvature of that accepted search
829 : !> \param position common line-search position used for comparison
830 : !> \param baseline_drop predicted conventional Mermin decrease
831 : !> \param response_drop predicted response Mermin decrease
832 : !> \param relative_gain response improvement relative to the larger predicted decrease
833 : !> \param preferred whether the response model predicts a resolved improvement
834 : ! **************************************************************************************************
835 32 : PURE SUBROUTINE ot_mermin_response_candidate_preferred( &
836 : baseline_slope, response_slope, response_curvature, accepted_slope, accepted_curvature, &
837 : position, baseline_drop, response_drop, relative_gain, preferred)
838 : REAL(KIND=dp), INTENT(IN) :: baseline_slope, response_slope, &
839 : response_curvature, accepted_slope, &
840 : accepted_curvature, position
841 : REAL(KIND=dp), INTENT(OUT) :: baseline_drop, response_drop, &
842 : relative_gain
843 : LOGICAL, INTENT(OUT) :: preferred
844 :
845 : REAL(KIND=dp) :: alpha, baseline_curvature, &
846 : curvature_ratio, resolution, scale
847 :
848 32 : baseline_drop = 0.0_dp
849 32 : response_drop = 0.0_dp
850 32 : relative_gain = 0.0_dp
851 32 : preferred = .FALSE.
852 : IF (.NOT. ieee_is_finite(baseline_slope) .OR. .NOT. ieee_is_finite(response_slope) .OR. &
853 : .NOT. ieee_is_finite(response_curvature) .OR. .NOT. ieee_is_finite(accepted_slope) .OR. &
854 32 : .NOT. ieee_is_finite(accepted_curvature) .OR. .NOT. ieee_is_finite(position)) RETURN
855 :
856 32 : alpha = ABS(position)
857 : scale = MAX(1.0_dp, ABS(baseline_slope), ABS(response_slope), &
858 32 : ABS(response_curvature), ABS(accepted_slope), ABS(accepted_curvature))
859 32 : resolution = 256.0_dp*EPSILON(1.0_dp)*scale
860 : IF (alpha <= SQRT(EPSILON(1.0_dp)) .OR. baseline_slope <= resolution .OR. &
861 32 : response_slope <= resolution .OR. accepted_slope <= resolution) RETURN
862 :
863 32 : curvature_ratio = accepted_curvature/accepted_slope
864 32 : baseline_curvature = curvature_ratio*baseline_slope
865 32 : baseline_drop = alpha*baseline_slope - 0.5_dp*alpha**2*baseline_curvature
866 32 : response_drop = alpha*response_slope - 0.5_dp*alpha**2*response_curvature
867 32 : scale = MAX(ABS(baseline_drop), ABS(response_drop), resolution)
868 32 : relative_gain = (response_drop - baseline_drop)/scale
869 44 : preferred = response_drop > resolution .AND. response_drop > baseline_drop + resolution
870 :
871 : END SUBROUTINE ot_mermin_response_candidate_preferred
872 :
873 : ! **************************************************************************************************
874 : !> \brief Select a sparse coupled-response probe from accepted-history evidence.
875 : !> \param available whether every local channel has a finite coupled response
876 : !> \param residual current conventional-preconditioned residual
877 : !> \param directions number of accepted directions seen by this response state
878 : !> \param shadow_good_samples passive accepted steps where response improved the prediction
879 : !> \param good_samples consecutive useful active response samples
880 : !> \param cooldown accepted directions remaining after a failed response sample
881 : !> \return whether the next direction should use the coupled response candidate
882 : ! **************************************************************************************************
883 23592 : PURE FUNCTION ot_mermin_response_probe( &
884 : available, residual, directions, shadow_good_samples, good_samples, cooldown) RESULT(probe)
885 : LOGICAL, INTENT(IN) :: available
886 : REAL(KIND=dp), INTENT(IN) :: residual
887 : INTEGER, INTENT(IN) :: directions, shadow_good_samples, &
888 : good_samples, cooldown
889 : LOGICAL :: probe
890 :
891 : INTEGER :: interval, phase
892 :
893 23592 : probe = .FALSE.
894 23592 : IF (.NOT. available .OR. .NOT. ieee_is_finite(residual)) RETURN
895 120 : IF (residual <= 0.0_dp .OR. residual > 2.0E-3_dp .OR. cooldown > 0) RETURN
896 120 : IF (shadow_good_samples < 3) RETURN
897 58 : interval = 8
898 58 : IF (good_samples >= 3) interval = 4
899 58 : IF (directions < interval) RETURN
900 36 : phase = MOD(directions, interval)
901 36 : probe = phase == 0 .OR. phase == 1
902 :
903 36 : END FUNCTION ot_mermin_response_probe
904 :
905 : ! **************************************************************************************************
906 : !> \brief Decide whether the next accepted state needs the dense finite Mermin response.
907 : !>
908 : !> A pending shadow is always completed. Otherwise the response is prepared while collecting
909 : !> the initial accepted-step calibration and at the state immediately preceding a sparse
910 : !> response-probe window. The relaxed residual bound accounts for using the preceding
911 : !> accepted residual before the current norm has been reduced across K-point groups.
912 : !> \param residual preceding accepted conventional-preconditioned residual
913 : !> \param directions number of accepted directions already seen
914 : !> \param shadow_good_samples useful passive accepted-step samples
915 : !> \param good_samples useful active response samples
916 : !> \param cooldown accepted directions remaining after a failed active sample
917 : !> \param shadow_pending whether the current endpoint must assess a prepared shadow
918 : !> \return whether to build the dense finite response at the current accepted endpoint
919 : ! **************************************************************************************************
920 396 : PURE FUNCTION ot_mermin_response_preparation_needed( &
921 : residual, directions, shadow_good_samples, good_samples, cooldown, shadow_pending) RESULT(prepare)
922 : REAL(KIND=dp), INTENT(IN) :: residual
923 : INTEGER, INTENT(IN) :: directions, shadow_good_samples, &
924 : good_samples, cooldown
925 : LOGICAL, INTENT(IN) :: shadow_pending
926 : LOGICAL :: prepare
927 :
928 : INTEGER :: interval, next_direction, phase
929 :
930 396 : prepare = shadow_pending
931 396 : IF (prepare) RETURN
932 252 : IF (.NOT. ieee_is_finite(residual) .OR. residual <= 0.0_dp .OR. residual > 4.0E-3_dp) RETURN
933 176 : IF (cooldown > 1) RETURN
934 160 : IF (shadow_good_samples < 3) THEN
935 396 : prepare = .TRUE.
936 : RETURN
937 : END IF
938 :
939 86 : interval = 8
940 86 : IF (good_samples >= 3) interval = 4
941 86 : next_direction = MIN(HUGE(directions) - 1, MAX(0, directions) + 1)
942 86 : phase = MOD(next_direction, interval)
943 86 : prepare = phase == interval - 1
944 :
945 86 : END FUNCTION ot_mermin_response_preparation_needed
946 :
947 : ! **************************************************************************************************
948 : !> \brief Decide whether a prepared conventional direction needs a shadow at its endpoint.
949 : !> \param residual current conventional-preconditioned residual
950 : !> \param directions number of accepted directions including the current one
951 : !> \param shadow_good_samples useful passive accepted-step samples
952 : !> \param good_samples useful active response samples
953 : !> \param cooldown accepted directions remaining after a failed active sample
954 : !> \return whether the next accepted endpoint must assess the current response shadow
955 : ! **************************************************************************************************
956 172 : PURE FUNCTION ot_mermin_response_shadow_followup( &
957 : residual, directions, shadow_good_samples, good_samples, cooldown) RESULT(followup)
958 : REAL(KIND=dp), INTENT(IN) :: residual
959 : INTEGER, INTENT(IN) :: directions, shadow_good_samples, &
960 : good_samples, cooldown
961 : LOGICAL :: followup
962 :
963 : INTEGER :: interval, phase
964 :
965 172 : followup = .FALSE.
966 172 : IF (.NOT. ieee_is_finite(residual) .OR. residual <= 0.0_dp .OR. residual > 4.0E-3_dp) RETURN
967 172 : IF (cooldown > 1) RETURN
968 172 : IF (shadow_good_samples < 3) THEN
969 172 : followup = .TRUE.
970 : RETURN
971 : END IF
972 :
973 60 : interval = 8
974 60 : IF (good_samples >= 3) interval = 4
975 60 : phase = MOD(MAX(0, directions), interval)
976 60 : followup = phase == interval - 1 .OR. phase == 0
977 :
978 : END FUNCTION ot_mermin_response_shadow_followup
979 :
980 : ! **************************************************************************************************
981 : !> \brief Assess and optionally select a calibrated finite Mermin response direction.
982 : !>
983 : !> The caller supplies its conventional descent direction in the product workspaces. A
984 : !> finite-response direction is considered only after accepted conventional steps have
985 : !> calibrated its quadratic model. This keeps the response independent of the line-search
986 : !> minimizer while preserving the physical residual as the convergence measure.
987 : !> \param qs_ot_env OT environments for all local spin/k-point channels
988 : !> \param para_env_inter_kp communicator between distributed k-point groups
989 : !> \param baseline_delta conventional preconditioned residual
990 : !> \param test_down physical gradient dotted into the conventional direction, updated on selection
991 : !> \param use_response_candidate whether the finite-response direction was selected
992 : ! **************************************************************************************************
993 23592 : SUBROUTINE ot_try_mermin_response_direction( &
994 : qs_ot_env, para_env_inter_kp, baseline_delta, test_down, use_response_candidate)
995 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
996 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
997 : REAL(KIND=dp), INTENT(IN) :: baseline_delta
998 : REAL(KIND=dp), INTENT(INOUT) :: test_down
999 : LOGICAL, INTENT(OUT) :: use_response_candidate
1000 :
1001 : INTEGER :: ispin, nspin
1002 : LOGICAL :: active_candidate_available, candidate_available, candidate_good, do_ener, do_ks, &
1003 : probe_response_candidate, shadow_good
1004 : REAL(KIND=dp) :: baseline_model_drop, measured_drop, predicted_drop, response_advantage, &
1005 : response_model_drop, response_quality, response_relative_gain, response_residual_ratio, &
1006 : response_test_down, tmp
1007 : TYPE(cp_logger_type), POINTER :: logger
1008 :
1009 23592 : nspin = SIZE(qs_ot_env)
1010 23592 : do_ks = qs_ot_env(1)%settings%ks
1011 23592 : do_ener = qs_ot_env(1)%settings%do_ener
1012 23592 : logger => cp_get_default_logger()
1013 23592 : shadow_good = .FALSE.
1014 23592 : candidate_good = .FALSE.
1015 :
1016 23592 : IF (qs_ot_env(1)%response_shadow_pending) THEN
1017 : CALL ot_mermin_response_compare( &
1018 : qs_ot_env(1)%response_reference_energy, qs_ot_env(1)%etotal, &
1019 : qs_ot_env(1)%response_reference_residual, baseline_delta, &
1020 : qs_ot_env(1)%response_predicted_slope, qs_ot_env(1)%response_shadow_curvature, &
1021 144 : qs_ot_env(1)%ds_min, response_advantage, response_residual_ratio, shadow_good)
1022 144 : IF (shadow_good) THEN
1023 : qs_ot_env(1)%response_shadow_good_samples = &
1024 : MIN(HUGE(qs_ot_env(1)%response_shadow_good_samples) - 1, &
1025 128 : qs_ot_env(1)%response_shadow_good_samples + 1)
1026 : ELSE
1027 : qs_ot_env(1)%response_shadow_good_samples = &
1028 16 : MAX(0, qs_ot_env(1)%response_shadow_good_samples - 1)
1029 : END IF
1030 144 : IF (logger%iter_info%print_level >= high_print_level .AND. logger%para_env%is_source()) THEN
1031 : WRITE (cp_logger_get_default_unit_nr(logger), &
1032 : '(A,1X,I5,2(1X,L1),4(1X,ES16.8))') &
1033 0 : " OT Mermin response shadow good hxc-dir advantage residual model shadow:", &
1034 0 : qs_ot_env(1)%response_shadow_good_samples, shadow_good, &
1035 0 : qs_ot_env(1)%response_hxc_direction_valid, response_advantage, &
1036 0 : response_residual_ratio, qs_ot_env(1)%response_model_curvature, &
1037 0 : qs_ot_env(1)%response_shadow_curvature
1038 : END IF
1039 144 : qs_ot_env(1)%response_shadow_pending = .FALSE.
1040 : END IF
1041 :
1042 23592 : IF (qs_ot_env(1)%response_candidate_pending) THEN
1043 : CALL ot_mermin_response_assess( &
1044 : qs_ot_env(1)%response_reference_energy, qs_ot_env(1)%etotal, &
1045 : qs_ot_env(1)%response_reference_residual, baseline_delta, &
1046 : qs_ot_env(1)%response_predicted_slope, qs_ot_env(1)%response_predicted_curvature, &
1047 : qs_ot_env(1)%ds_min, qs_ot_env(1)%settings%ds_min, predicted_drop, measured_drop, &
1048 16 : response_quality, response_residual_ratio, candidate_good)
1049 16 : IF (candidate_good) THEN
1050 : qs_ot_env(1)%response_candidate_good_samples = &
1051 : MIN(HUGE(qs_ot_env(1)%response_candidate_good_samples) - 1, &
1052 12 : qs_ot_env(1)%response_candidate_good_samples + 1)
1053 : ELSE
1054 : qs_ot_env(1)%response_candidate_good_samples = &
1055 4 : MAX(0, qs_ot_env(1)%response_candidate_good_samples - 2)
1056 4 : qs_ot_env(1)%response_candidate_cooldown = 6
1057 : END IF
1058 16 : IF (logger%iter_info%print_level >= high_print_level .AND. logger%para_env%is_source()) THEN
1059 : WRITE (cp_logger_get_default_unit_nr(logger), '(A,1X,I5,1X,L1,4(1X,ES16.8))') &
1060 0 : " OT Mermin response candidate good quality residual predicted measured:", &
1061 0 : qs_ot_env(1)%response_candidate_good_samples, candidate_good, response_quality, &
1062 0 : response_residual_ratio, predicted_drop, measured_drop
1063 : END IF
1064 16 : qs_ot_env(1)%response_candidate_pending = .FALSE.
1065 : END IF
1066 :
1067 : qs_ot_env(1)%response_candidate_directions = &
1068 : MIN(HUGE(qs_ot_env(1)%response_candidate_directions) - 1, &
1069 23592 : qs_ot_env(1)%response_candidate_directions + 1)
1070 23592 : IF (qs_ot_env(1)%response_candidate_cooldown > 0) THEN
1071 24 : qs_ot_env(1)%response_candidate_cooldown = qs_ot_env(1)%response_candidate_cooldown - 1
1072 : END IF
1073 :
1074 : candidate_available = qs_ot_env(1)%settings%occupation_preconditioner .AND. &
1075 23592 : qs_ot_env(1)%settings%do_rotation .AND. do_ener
1076 52838 : DO ispin = 1, nspin
1077 : candidate_available = candidate_available .AND. &
1078 : qs_ot_env(ispin)%rotation_response_valid .AND. &
1079 : ASSOCIATED(qs_ot_env(ispin)%matrix_response_gx) .AND. &
1080 : ASSOCIATED(qs_ot_env(ispin)%matrix_response_gx_im) .AND. &
1081 : ASSOCIATED(qs_ot_env(ispin)%rot_mat_response_gx) .AND. &
1082 : ASSOCIATED(qs_ot_env(ispin)%rot_mat_response_gx_im) .AND. &
1083 81814 : ASSOCIATED(qs_ot_env(ispin)%ener_response_gx)
1084 : END DO
1085 : active_candidate_available = candidate_available .AND. shadow_good .AND. &
1086 23592 : qs_ot_env(1)%response_hxc_direction_valid
1087 : probe_response_candidate = ot_mermin_response_probe( &
1088 : active_candidate_available, baseline_delta, &
1089 : qs_ot_env(1)%response_candidate_directions, &
1090 : qs_ot_env(1)%response_shadow_good_samples, &
1091 : qs_ot_env(1)%response_candidate_good_samples, &
1092 23592 : qs_ot_env(1)%response_candidate_cooldown)
1093 23592 : use_response_candidate = .FALSE.
1094 :
1095 23592 : IF (probe_response_candidate) THEN
1096 32 : response_test_down = 0.0_dp
1097 32 : baseline_model_drop = 0.0_dp
1098 32 : response_model_drop = 0.0_dp
1099 32 : response_relative_gain = 0.0_dp
1100 32 : IF (do_ks) THEN
1101 76 : DO ispin = 1, nspin
1102 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
1103 44 : qs_ot_env(ispin)%matrix_response_gx, tmp)
1104 44 : response_test_down = response_test_down - tmp
1105 44 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1106 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
1107 44 : qs_ot_env(ispin)%matrix_response_gx_im, tmp)
1108 44 : response_test_down = response_test_down - tmp
1109 : END IF
1110 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
1111 44 : qs_ot_env(ispin)%rot_mat_response_gx, tmp)
1112 44 : response_test_down = response_test_down - 0.5_dp*tmp
1113 120 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1114 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
1115 44 : qs_ot_env(ispin)%rot_mat_response_gx_im, tmp)
1116 44 : response_test_down = response_test_down - 0.5_dp*tmp
1117 : END IF
1118 : END DO
1119 : END IF
1120 32 : IF (do_ener) THEN
1121 76 : DO ispin = 1, nspin
1122 : response_test_down = response_test_down - &
1123 : DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
1124 292 : qs_ot_env(ispin)%ener_response_gx)
1125 : END DO
1126 : END IF
1127 32 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, response_test_down)
1128 32 : IF (response_test_down < 0.0_dp) THEN
1129 : CALL ot_mermin_response_candidate_preferred( &
1130 : -test_down, -response_test_down, qs_ot_env(1)%response_model_curvature, &
1131 : qs_ot_env(1)%response_predicted_slope, qs_ot_env(1)%response_shadow_curvature, &
1132 : qs_ot_env(1)%ds_min, baseline_model_drop, response_model_drop, &
1133 32 : response_relative_gain, use_response_candidate)
1134 : END IF
1135 32 : IF (logger%iter_info%print_level >= high_print_level .AND. logger%para_env%is_source()) THEN
1136 : WRITE (cp_logger_get_default_unit_nr(logger), '(A,4(1X,I5),5(1X,ES16.8),2(1X,L1))') &
1137 : " OT Mermin response candidate comparison directions shadow active cooldown residual "// &
1138 0 : "curvature baseline response gain probe use:", &
1139 0 : qs_ot_env(1)%response_candidate_directions, &
1140 0 : qs_ot_env(1)%response_shadow_good_samples, &
1141 0 : qs_ot_env(1)%response_candidate_good_samples, &
1142 0 : qs_ot_env(1)%response_candidate_cooldown, baseline_delta, &
1143 0 : qs_ot_env(1)%response_model_curvature, baseline_model_drop, response_model_drop, &
1144 0 : response_relative_gain, probe_response_candidate, use_response_candidate
1145 : END IF
1146 32 : IF (use_response_candidate) THEN
1147 48 : DO ispin = 1, nspin
1148 28 : IF (do_ks) THEN
1149 28 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_response_gx)
1150 28 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
1151 28 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1152 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx_im, &
1153 28 : qs_ot_env(ispin)%matrix_response_gx_im)
1154 28 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx_im, -1.0_dp)
1155 : END IF
1156 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, &
1157 28 : qs_ot_env(ispin)%rot_mat_response_gx)
1158 28 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
1159 28 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1160 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, &
1161 28 : qs_ot_env(ispin)%rot_mat_response_gx_im)
1162 28 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx_im, -1.0_dp)
1163 : END IF
1164 : END IF
1165 184 : IF (do_ener) qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_response_gx
1166 : END DO
1167 20 : test_down = response_test_down
1168 : END IF
1169 : END IF
1170 :
1171 23592 : qs_ot_env(1)%response_reference_energy = qs_ot_env(1)%etotal
1172 23592 : qs_ot_env(1)%response_reference_residual = baseline_delta
1173 23592 : IF (test_down < 0.0_dp) THEN
1174 22920 : qs_ot_env(1)%response_predicted_slope = -test_down
1175 : ELSE
1176 672 : qs_ot_env(1)%response_predicted_slope = qs_ot_env(1)%gnorm
1177 : END IF
1178 23592 : IF (use_response_candidate) THEN
1179 20 : qs_ot_env(1)%response_candidate_pending = .TRUE.
1180 20 : qs_ot_env(1)%response_shadow_pending = .FALSE.
1181 20 : qs_ot_env(1)%response_predicted_curvature = qs_ot_env(1)%response_model_curvature
1182 : ELSE
1183 23572 : qs_ot_env(1)%response_candidate_pending = .FALSE.
1184 : qs_ot_env(1)%response_shadow_pending = candidate_available .AND. &
1185 : ot_mermin_response_shadow_followup( &
1186 : baseline_delta, &
1187 : qs_ot_env(1)%response_candidate_directions, &
1188 : qs_ot_env(1)%response_shadow_good_samples, &
1189 : qs_ot_env(1)%response_candidate_good_samples, &
1190 46998 : qs_ot_env(1)%response_candidate_cooldown)
1191 : END IF
1192 :
1193 23592 : END SUBROUTINE ot_try_mermin_response_direction
1194 :
1195 : !
1196 : ! checks if the current point is a good point for finding a new direction
1197 : ! or if we should improve the line_search, if it is used
1198 : !
1199 : ! **************************************************************************************************
1200 : !> \brief ...
1201 : !> \param qs_ot_env ...
1202 : !> \return ...
1203 : ! **************************************************************************************************
1204 45540 : FUNCTION current_point_is_fine(qs_ot_env) RESULT(res)
1205 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1206 : LOGICAL :: res
1207 :
1208 45540 : res = .FALSE.
1209 :
1210 : ! only if we have a gradient it can be fine
1211 45540 : IF (.NOT. qs_ot_env(1)%energy_only) THEN
1212 :
1213 : ! we have not yet started with the line search
1214 23642 : IF (qs_ot_env(1)%line_search_count == 0) THEN
1215 45540 : res = .TRUE.
1216 : RETURN
1217 : END IF
1218 :
1219 20396 : IF (qs_ot_env(1)%line_search_might_be_done) THEN
1220 : ! here we put the more complicated logic later
1221 20396 : res = .TRUE.
1222 20396 : RETURN
1223 : END IF
1224 :
1225 : END IF
1226 :
1227 : END FUNCTION current_point_is_fine
1228 :
1229 : !
1230 : ! performs various kinds of line searches
1231 : !
1232 : ! **************************************************************************************************
1233 : !> \brief ...
1234 : !> \param qs_ot_env ...
1235 : ! **************************************************************************************************
1236 45540 : SUBROUTINE do_line_search(qs_ot_env)
1237 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1238 :
1239 45968 : SELECT CASE (qs_ot_env(1)%settings%line_search_method)
1240 : CASE ("GOLD")
1241 428 : CALL do_line_search_gold(qs_ot_env)
1242 : CASE ("3PNT")
1243 1656 : CALL do_line_search_3pnt(qs_ot_env)
1244 : CASE ("2PNT")
1245 43366 : IF (use_three_point_mermin_search(qs_ot_env)) THEN
1246 1136 : CALL do_line_search_3pnt(qs_ot_env)
1247 : ELSE
1248 42230 : CALL do_line_search_2pnt(qs_ot_env)
1249 : END IF
1250 : CASE ("ADPT")
1251 80 : CALL do_line_search_adapt(qs_ot_env)
1252 : CASE ("NONE")
1253 10 : CALL do_line_search_none(qs_ot_env)
1254 : CASE DEFAULT
1255 45540 : CPABORT("NYI")
1256 : END SELECT
1257 45540 : END SUBROUTINE do_line_search
1258 :
1259 : ! **************************************************************************************************
1260 : !> \brief Add an energy-only guard point to strongly preconditioned complex Mermin steps.
1261 : !> \param qs_ot_env OT environments
1262 : !> \return Whether to use the guarded three-point search
1263 : ! **************************************************************************************************
1264 44270 : FUNCTION use_three_point_mermin_search(qs_ot_env) RESULT(res)
1265 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1266 : LOGICAL :: res
1267 :
1268 : INTEGER :: ispin
1269 :
1270 44270 : res = .FALSE.
1271 44270 : IF (.NOT. qs_ot_env(1)%settings%occupation_preconditioner) RETURN
1272 :
1273 1500 : DO ispin = 1, SIZE(qs_ot_env)
1274 1500 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1275 44270 : res = .TRUE.
1276 : RETURN
1277 : END IF
1278 : END DO
1279 : END FUNCTION use_three_point_mermin_search
1280 :
1281 : ! **************************************************************************************************
1282 : !> \brief moves x adding the right amount (ds) of the gradient or search direction
1283 : !> \param ds ...
1284 : !> \param qs_ot_env ...
1285 : !> \par History
1286 : !> 08.2004 created [ Joost VandeVondele ] copied here from a larger number of subroutines
1287 : ! **************************************************************************************************
1288 45540 : SUBROUTINE take_step(ds, qs_ot_env)
1289 : REAL(KIND=dp) :: ds
1290 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1291 :
1292 : CHARACTER(len=*), PARAMETER :: routineN = 'take_step'
1293 :
1294 : INTEGER :: handle, ispin, nspin
1295 : LOGICAL :: do_ener, do_ks
1296 :
1297 45540 : CALL timeset(routineN, handle)
1298 :
1299 45540 : nspin = SIZE(qs_ot_env)
1300 :
1301 45540 : do_ks = qs_ot_env(1)%settings%ks
1302 45540 : do_ener = qs_ot_env(1)%settings%do_ener
1303 :
1304 : ! now update x to take into account this new step
1305 : ! either dx or -gx is the direction to use
1306 45540 : IF (qs_ot_env(1)%use_dx) THEN
1307 45452 : IF (do_ks) THEN
1308 101460 : DO ispin = 1, nspin
1309 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_dx, &
1310 56008 : alpha_scalar=1.0_dp, beta_scalar=ds)
1311 56008 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1312 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x_im, qs_ot_env(ispin)%matrix_dx_im, &
1313 3603 : alpha_scalar=1.0_dp, beta_scalar=ds)
1314 : END IF
1315 101460 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1316 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, qs_ot_env(ispin)%rot_mat_dx, &
1317 4963 : alpha_scalar=1.0_dp, beta_scalar=ds)
1318 4963 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1319 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x_im, qs_ot_env(ispin)%rot_mat_dx_im, &
1320 2369 : alpha_scalar=1.0_dp, beta_scalar=ds)
1321 : END IF
1322 : END IF
1323 : END DO
1324 : END IF
1325 45452 : IF (do_ener) THEN
1326 2982 : DO ispin = 1, nspin
1327 14510 : qs_ot_env(ispin)%ener_x = qs_ot_env(ispin)%ener_x + ds*qs_ot_env(ispin)%ener_dx
1328 : END DO
1329 : END IF
1330 : ELSE
1331 88 : IF (do_ks) THEN
1332 176 : DO ispin = 1, nspin
1333 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_gx, &
1334 88 : alpha_scalar=1.0_dp, beta_scalar=-ds)
1335 88 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1336 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x_im, qs_ot_env(ispin)%matrix_gx_im, &
1337 0 : alpha_scalar=1.0_dp, beta_scalar=-ds)
1338 : END IF
1339 176 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1340 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, qs_ot_env(ispin)%rot_mat_gx, &
1341 0 : alpha_scalar=1.0_dp, beta_scalar=-ds)
1342 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1343 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x_im, qs_ot_env(ispin)%rot_mat_gx_im, &
1344 0 : alpha_scalar=1.0_dp, beta_scalar=-ds)
1345 : END IF
1346 : END IF
1347 : END DO
1348 : END IF
1349 88 : IF (do_ener) THEN
1350 0 : DO ispin = 1, nspin
1351 0 : qs_ot_env(ispin)%ener_x = qs_ot_env(ispin)%ener_x - ds*qs_ot_env(ispin)%ener_gx
1352 : END DO
1353 : END IF
1354 : END IF
1355 45540 : CALL timestop(handle)
1356 45540 : END SUBROUTINE take_step
1357 :
1358 : ! implements a golden ratio search as a robust way of minimizing
1359 : ! **************************************************************************************************
1360 : !> \brief ...
1361 : !> \param qs_ot_env ...
1362 : ! **************************************************************************************************
1363 428 : SUBROUTINE do_line_search_gold(qs_ot_env)
1364 :
1365 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1366 :
1367 : CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_gold'
1368 : REAL(KIND=dp), PARAMETER :: gold_sec = 0.3819_dp
1369 :
1370 : INTEGER :: count, handle
1371 : REAL(KIND=dp) :: ds
1372 :
1373 428 : CALL timeset(routineN, handle)
1374 :
1375 428 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
1376 428 : count = qs_ot_env(1)%line_search_count
1377 428 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
1378 428 : qs_ot_env(1)%energy_only = .TRUE.
1379 :
1380 428 : IF (count + 1 > SIZE(qs_ot_env(1)%OT_pos)) THEN
1381 : ! should not happen, we pass with a warning first
1382 : ! you can increase the size of OT_pos and the like in qs_ot_env
1383 0 : CPABORT("MAX ITER EXCEEDED : FATAL")
1384 : END IF
1385 :
1386 428 : IF (qs_ot_env(1)%line_search_count == 1) THEN
1387 54 : qs_ot_env(1)%line_search_left = 1
1388 54 : qs_ot_env(1)%line_search_right = 0
1389 54 : qs_ot_env(1)%line_search_mid = 1
1390 54 : qs_ot_env(1)%ot_pos(1) = 0.0_dp
1391 54 : qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
1392 54 : qs_ot_env(1)%ot_pos(2) = qs_ot_env(1)%ds_min/gold_sec
1393 : ELSE
1394 374 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
1395 : ! it's essentially a book keeping game.
1396 : ! keep left on the left, keep (bring) right on the right
1397 : ! and mid in between these two
1398 374 : IF (qs_ot_env(1)%line_search_right == 0) THEN ! we do not yet have the right bracket
1399 64 : IF (qs_ot_env(1)%ot_energy(count - 1) < qs_ot_env(1)%ot_energy(count)) THEN
1400 42 : qs_ot_env(1)%line_search_right = count
1401 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) + &
1402 : (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) - &
1403 42 : qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))*gold_sec
1404 : ELSE
1405 22 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
1406 22 : qs_ot_env(1)%line_search_mid = count
1407 22 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ot_pos(count)/gold_sec ! expand
1408 : END IF
1409 : ELSE
1410 : ! first determine where we are and construct the new triplet
1411 310 : IF (qs_ot_env(1)%ot_pos(count) < qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
1412 126 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
1413 42 : qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
1414 42 : qs_ot_env(1)%line_search_mid = count
1415 : ELSE
1416 84 : qs_ot_env(1)%line_search_left = count
1417 : END IF
1418 : ELSE
1419 184 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
1420 76 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
1421 76 : qs_ot_env(1)%line_search_mid = count
1422 : ELSE
1423 108 : qs_ot_env(1)%line_search_right = count
1424 : END IF
1425 : END IF
1426 : ! now find the new point in the largest section
1427 310 : IF ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
1428 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) > &
1429 : (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
1430 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left))) THEN
1431 : qs_ot_env(1)%ot_pos(count + 1) = &
1432 : qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) + &
1433 : gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
1434 162 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))
1435 : ELSE
1436 : qs_ot_env(1)%ot_pos(count + 1) = &
1437 : qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left) + &
1438 : gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
1439 148 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left))
1440 : END IF
1441 : ! check for termination
1442 : IF (((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
1443 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) < &
1444 310 : qs_ot_env(1)%ds_min*qs_ot_env(1)%settings%gold_target) .AND. &
1445 : ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
1446 : - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left)) < &
1447 : qs_ot_env(1)%ds_min*qs_ot_env(1)%settings%gold_target)) THEN
1448 42 : qs_ot_env(1)%energy_only = .FALSE.
1449 42 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
1450 : END IF
1451 : END IF
1452 : END IF
1453 428 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
1454 428 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
1455 :
1456 428 : CALL take_step(ds, qs_ot_env)
1457 :
1458 428 : CALL timestop(handle)
1459 :
1460 428 : END SUBROUTINE do_line_search_gold
1461 :
1462 : ! **************************************************************************************************
1463 : !> \brief ...
1464 : !> \param qs_ot_env ...
1465 : ! **************************************************************************************************
1466 2792 : SUBROUTINE do_line_search_3pnt(qs_ot_env)
1467 :
1468 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1469 :
1470 : CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_3pnt'
1471 :
1472 : INTEGER :: best_index, count, handle, i
1473 : REAL(KIND=dp) :: best_energy, denom, ds, fa, fb, fc, nom, &
1474 : pos, tol, val, xa, xb, xc
1475 :
1476 2792 : CALL timeset(routineN, handle)
1477 :
1478 2792 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
1479 2792 : qs_ot_env(1)%energy_only = .TRUE.
1480 :
1481 : ! a three point interpolation based on the energy
1482 2792 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
1483 2792 : count = qs_ot_env(1)%line_search_count
1484 2792 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
1485 984 : SELECT CASE (count)
1486 : CASE (1)
1487 984 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
1488 984 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*0.8_dp
1489 : CASE (2)
1490 904 : IF (qs_ot_env(1)%OT_energy(count) > qs_ot_env(1)%OT_energy(count - 1)) THEN
1491 88 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*0.5_dp
1492 : ELSE
1493 816 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*1.4_dp
1494 : END IF
1495 : CASE (3)
1496 904 : xa = qs_ot_env(1)%OT_pos(1)
1497 904 : xb = qs_ot_env(1)%OT_pos(2)
1498 904 : xc = qs_ot_env(1)%OT_pos(3)
1499 904 : fa = qs_ot_env(1)%OT_energy(1)
1500 904 : fb = qs_ot_env(1)%OT_energy(2)
1501 904 : fc = qs_ot_env(1)%OT_energy(3)
1502 904 : nom = (xb - xa)**2*(fb - fc) - (xb - xc)**2*(fb - fa)
1503 904 : denom = (xb - xa)*(fb - fc) - (xb - xc)*(fb - fa)
1504 904 : IF (ABS(denom) <= 1.0E-18_dp*MAX(ABS(fb - fc), ABS(fb - fa))) THEN
1505 : pos = xb
1506 : ELSE
1507 904 : pos = xb - 0.5_dp*nom/denom ! position of the stationary point
1508 : END IF
1509 : val = (pos - xa)*(pos - xb)*fc/((xc - xa)*(xc - xb)) + &
1510 : (pos - xb)*(pos - xc)*fa/((xa - xb)*(xa - xc)) + &
1511 904 : (pos - xc)*(pos - xa)*fb/((xb - xc)*(xb - xa))
1512 904 : best_index = 1
1513 2712 : DO i = 2, count
1514 2712 : IF (qs_ot_env(1)%OT_energy(i) < qs_ot_env(1)%OT_energy(best_index)) best_index = i
1515 : END DO
1516 904 : best_energy = qs_ot_env(1)%OT_energy(best_index)
1517 904 : tol = 10.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(best_energy))
1518 904 : IF (use_three_point_mermin_search(qs_ot_env) .AND. best_index /= 1 .AND. &
1519 : val >= best_energy - tol) THEN
1520 6 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%OT_pos(best_index)
1521 898 : ELSE IF (val < fa .AND. val <= fb .AND. val <= fc) THEN ! OK, we go to a minimum
1522 : ! we take a guard against too large steps
1523 : qs_ot_env(1)%OT_pos(count + 1) = MAX(MAXVAL(qs_ot_env(1)%OT_pos(1:3))*0.01_dp, &
1524 3584 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:3))*4.0_dp))
1525 : ELSE ! just take an extended step
1526 8 : qs_ot_env(1)%OT_pos(count + 1) = MAXVAL(qs_ot_env(1)%OT_pos(1:3))*2.0_dp
1527 : END IF
1528 904 : qs_ot_env(1)%energy_only = .FALSE.
1529 904 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
1530 : CASE DEFAULT
1531 2792 : CPABORT("NYI")
1532 : END SELECT
1533 2792 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
1534 2792 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
1535 :
1536 2792 : CALL take_step(ds, qs_ot_env)
1537 :
1538 2792 : CALL timestop(handle)
1539 :
1540 2792 : END SUBROUTINE do_line_search_3pnt
1541 :
1542 : ! **************************************************************************************************
1543 : !> \brief ...
1544 : !> \param qs_ot_env ...
1545 : ! **************************************************************************************************
1546 42230 : SUBROUTINE do_line_search_2pnt(qs_ot_env)
1547 :
1548 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1549 :
1550 : CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_2pnt'
1551 :
1552 : INTEGER :: count, handle
1553 : REAL(KIND=dp) :: a, b, c, ds, pos, val, x0, x1
1554 :
1555 42230 : CALL timeset(routineN, handle)
1556 :
1557 42230 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
1558 42230 : qs_ot_env(1)%energy_only = .TRUE.
1559 :
1560 : ! a three point interpolation based on the energy
1561 42230 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
1562 42230 : count = qs_ot_env(1)%line_search_count
1563 42230 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
1564 22574 : SELECT CASE (count)
1565 : CASE (1)
1566 22574 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
1567 22574 : qs_ot_env(1)%ot_grad(count) = qs_ot_env(1)%gradient
1568 22574 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*1.0_dp
1569 : CASE (2)
1570 19656 : x0 = 0.0_dp
1571 19656 : c = qs_ot_env(1)%ot_energy(1)
1572 19656 : b = qs_ot_env(1)%ot_grad(1)
1573 19656 : x1 = qs_ot_env(1)%ot_pos(2)
1574 19656 : a = (qs_ot_env(1)%ot_energy(2) - b*x1 - c)/(x1**2)
1575 19656 : IF (a <= 0.0_dp) a = 1.0E-15_dp
1576 19656 : pos = -b/(2.0_dp*a)
1577 19656 : val = a*pos**2 + b*pos + c
1578 19656 : qs_ot_env(1)%energy_only = .FALSE.
1579 19656 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
1580 19656 : IF (val < qs_ot_env(1)%ot_energy(1) .AND. val <= qs_ot_env(1)%ot_energy(2)) THEN
1581 : ! we go to a minimum, but ...
1582 : ! we take a guard against too large steps
1583 : qs_ot_env(1)%OT_pos(count + 1) = MAX(MAXVAL(qs_ot_env(1)%OT_pos(1:2))*0.01_dp, &
1584 58758 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:2))*4.0_dp))
1585 : ELSE ! just take an extended step
1586 210 : qs_ot_env(1)%OT_pos(count + 1) = MAXVAL(qs_ot_env(1)%OT_pos(1:2))*2.0_dp
1587 : END IF
1588 : CASE DEFAULT
1589 42230 : CPABORT("NYI")
1590 : END SELECT
1591 42230 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
1592 42230 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
1593 :
1594 42230 : CALL take_step(ds, qs_ot_env)
1595 :
1596 42230 : CALL timestop(handle)
1597 :
1598 42230 : END SUBROUTINE do_line_search_2pnt
1599 :
1600 : ! **************************************************************************************************
1601 : !> \brief ...
1602 : !> \param qs_ot_env ...
1603 : ! **************************************************************************************************
1604 80 : SUBROUTINE do_line_search_adapt(qs_ot_env)
1605 :
1606 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1607 :
1608 : CHARACTER(len=*), PARAMETER :: routineN = 'do_line_search_adapt'
1609 : REAL(KIND=dp), PARAMETER :: grow_factor = 2.0_dp, &
1610 : shrink_factor = 0.5_dp
1611 :
1612 : INTEGER :: count, handle, il, im, ir
1613 : REAL(KIND=dp) :: a, b, c, denom, ds, el, em, er, &
1614 : step_size, xl, xm, xr
1615 :
1616 80 : CALL timeset(routineN, handle)
1617 :
1618 80 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
1619 80 : count = qs_ot_env(1)%line_search_count
1620 80 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
1621 80 : qs_ot_env(1)%energy_only = .TRUE.
1622 :
1623 80 : IF (count + 1 > SIZE(qs_ot_env(1)%OT_pos)) THEN
1624 : ! should not happen, we pass with a warning first
1625 : ! you can increase the size of OT_pos and the like in qs_ot_env
1626 0 : CPABORT("MAX ITER EXCEEDED : FATAL")
1627 : END IF
1628 :
1629 : ! Perform an adaptive linesearch
1630 80 : IF (qs_ot_env(1)%line_search_count == 1) THEN
1631 20 : qs_ot_env(1)%line_search_left = 1
1632 20 : qs_ot_env(1)%line_search_right = 0
1633 20 : qs_ot_env(1)%line_search_mid = 1
1634 20 : qs_ot_env(1)%ot_Pos(1) = 0.0_dp
1635 20 : qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
1636 20 : qs_ot_env(1)%ot_Pos(2) = qs_ot_env(1)%ds_min*grow_factor
1637 : ELSE
1638 60 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
1639 : ! it's essentially a book keeping game.
1640 : ! keep left on the left, keep (bring) right on the right
1641 : ! and mid in between these two
1642 60 : IF (qs_ot_env(1)%line_search_right == 0) THEN ! we do not yet have the right bracket
1643 28 : IF (qs_ot_env(1)%ot_energy(count - 1) < qs_ot_env(1)%ot_energy(count)) THEN
1644 16 : qs_ot_env(1)%line_search_right = count
1645 : qs_ot_env(1)%ot_Pos(count + 1) = qs_ot_env(1)%ot_Pos(qs_ot_env(1)%line_search_mid) + &
1646 : (qs_ot_env(1)%ot_Pos(qs_ot_env(1)%line_search_right) - &
1647 16 : qs_ot_env(1)%ot_Pos(qs_ot_env(1)%line_search_mid))*shrink_factor
1648 : ELSE
1649 : ! expand further
1650 12 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
1651 12 : qs_ot_env(1)%line_search_mid = count
1652 12 : qs_ot_env(1)%ot_Pos(count + 1) = qs_ot_env(1)%ot_Pos(count)*grow_factor
1653 : END IF
1654 : ELSE
1655 : ! first determine where we are and construct the new triplet
1656 32 : IF (qs_ot_env(1)%ot_pos(count) < qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
1657 0 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
1658 0 : qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
1659 0 : qs_ot_env(1)%line_search_mid = count
1660 : ELSE
1661 0 : qs_ot_env(1)%line_search_left = count
1662 : END IF
1663 : ELSE
1664 32 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
1665 12 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
1666 12 : qs_ot_env(1)%line_search_mid = count
1667 : ELSE
1668 20 : qs_ot_env(1)%line_search_right = count
1669 : END IF
1670 : END IF
1671 32 : il = qs_ot_env(1)%line_search_left
1672 32 : im = qs_ot_env(1)%line_search_mid
1673 32 : ir = qs_ot_env(1)%line_search_right
1674 32 : xl = qs_ot_env(1)%OT_pos(il)
1675 32 : xm = qs_ot_env(1)%OT_pos(im)
1676 32 : xr = qs_ot_env(1)%OT_pos(ir)
1677 32 : el = qs_ot_env(1)%ot_energy(il)
1678 32 : em = qs_ot_env(1)%ot_energy(im)
1679 32 : er = qs_ot_env(1)%ot_energy(ir)
1680 32 : IF (em < el) THEN
1681 16 : IF (er < em) THEN
1682 : !extend search
1683 0 : qs_ot_env(1)%ot_Pos(count + 1) = qs_ot_env(1)%ot_Pos(ir)*grow_factor
1684 : ELSE
1685 : ! Cramer's rule
1686 16 : denom = (xl - xm)*(xl - xr)*(xm - xr)
1687 16 : a = (xr*(em - el) + xm*(el - er) + xl*(er - em))/denom
1688 16 : b = (xr**2*(el - em) + xm**2*(er - el) + xl**2*(em - er))/denom
1689 16 : c = (xm*xr*(xm - xr)*el + xr*xl*(xr - xl)*em + xr*xm*(xr - xm)*er)/denom
1690 :
1691 16 : IF (ABS(a) /= 0.0_dp) THEN
1692 16 : step_size = -b/(2.0_dp*a)
1693 : ELSE
1694 : step_size = 0.0_dp
1695 : END IF
1696 16 : CPASSERT(step_size >= 0.0_dp)
1697 16 : qs_ot_env(1)%ot_Pos(count + 1) = step_size
1698 16 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
1699 16 : qs_ot_env(1)%energy_only = .FALSE.
1700 : END IF
1701 : ELSE
1702 : ! contract search
1703 : qs_ot_env(1)%ot_Pos(count + 1) = qs_ot_env(1)%ot_Pos(im) + &
1704 16 : (qs_ot_env(1)%ot_Pos(ir) - qs_ot_env(1)%ot_Pos(im))*shrink_factor
1705 : END IF
1706 :
1707 : END IF
1708 : END IF
1709 80 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
1710 80 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
1711 :
1712 80 : CALL take_step(ds, qs_ot_env)
1713 :
1714 80 : CALL timestop(handle)
1715 :
1716 80 : END SUBROUTINE do_line_search_adapt
1717 :
1718 : ! **************************************************************************************************
1719 : !> \brief ...
1720 : !> \param qs_ot_env ...
1721 : ! **************************************************************************************************
1722 10 : SUBROUTINE do_line_search_none(qs_ot_env)
1723 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1724 :
1725 10 : CALL take_step(qs_ot_env(1)%ds_min, qs_ot_env)
1726 :
1727 10 : END SUBROUTINE do_line_search_none
1728 :
1729 : !
1730 : ! creates a new SD direction, using the preconditioner if associated
1731 : ! also updates the gradient for line search
1732 : !
1733 :
1734 : ! **************************************************************************************************
1735 : !> \brief ...
1736 : !> \param qs_ot_env ...
1737 : !> \param para_env_inter_kp communicator between distributed k-point groups
1738 : ! **************************************************************************************************
1739 50 : SUBROUTINE ot_new_sd_direction(qs_ot_env, para_env_inter_kp)
1740 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1741 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
1742 :
1743 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_new_sd_direction'
1744 :
1745 : INTEGER :: handle, ispin, nspin
1746 : LOGICAL :: do_ener, do_ks, use_occupation_response
1747 : REAL(KIND=dp) :: nvariables, tmp
1748 : TYPE(cp_logger_type), POINTER :: logger
1749 : TYPE(dbcsr_type), POINTER :: orbital_residual, orbital_residual_im
1750 236 : TYPE(dbcsr_p_type) :: matrix_out(2, SIZE(qs_ot_env))
1751 :
1752 50 : CALL timeset(routineN, handle)
1753 :
1754 50 : nspin = SIZE(qs_ot_env)
1755 50 : do_ks = qs_ot_env(1)%settings%ks
1756 50 : do_ener = qs_ot_env(1)%settings%do_ener
1757 50 : use_occupation_response = qs_ot_env(1)%settings%occupation_preconditioner
1758 :
1759 50 : IF (do_ks .AND. qs_ot_env(1)%use_dx) THEN
1760 88 : DO ispin = 1, SIZE(qs_ot_env)
1761 50 : matrix_out(1, ispin)%matrix => qs_ot_env(ispin)%matrix_dx
1762 88 : matrix_out(2, ispin)%matrix => qs_ot_env(ispin)%matrix_dx_im
1763 : END DO
1764 : END IF
1765 :
1766 : ! Keep the physical gradients until the complete orbital/rotation/energy product is known.
1767 0 : DO
1768 50 : qs_ot_env(1)%gnorm = 0.0_dp
1769 50 : IF (do_ks) THEN
1770 50 : IF (qs_ot_env(1)%use_dx) CALL ot_precondition_gradients(qs_ot_env, matrix_out, use_occupation_response)
1771 112 : DO ispin = 1, nspin
1772 62 : IF (use_occupation_response) THEN
1773 24 : orbital_residual => qs_ot_env(ispin)%matrix_preconditioned_gx
1774 24 : orbital_residual_im => qs_ot_env(ispin)%matrix_preconditioned_gx_im
1775 : ELSE
1776 38 : orbital_residual => qs_ot_env(ispin)%matrix_gx
1777 38 : orbital_residual_im => qs_ot_env(ispin)%matrix_gx_im
1778 : END IF
1779 62 : IF (qs_ot_env(1)%use_dx) THEN
1780 50 : orbital_residual => qs_ot_env(ispin)%matrix_dx
1781 50 : orbital_residual_im => qs_ot_env(ispin)%matrix_dx_im
1782 : END IF
1783 62 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, orbital_residual, tmp)
1784 62 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1785 112 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1786 24 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, orbital_residual_im, tmp)
1787 24 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1788 : END IF
1789 : END DO
1790 50 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1791 36 : DO ispin = 1, nspin
1792 24 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
1793 : ! The antisymmetric rotation stores each variable twice.
1794 24 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
1795 36 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1796 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
1797 24 : qs_ot_env(ispin)%rot_mat_gx_im, tmp)
1798 24 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
1799 : END IF
1800 : END DO
1801 : END IF
1802 : END IF
1803 50 : IF (do_ener) THEN
1804 36 : DO ispin = 1, nspin
1805 24 : IF (use_occupation_response) THEN
1806 72 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_preconditioned_gx)
1807 : ELSE
1808 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
1809 : END IF
1810 36 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1811 : END DO
1812 : END IF
1813 50 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%gnorm)
1814 :
1815 : ! Match CG and history residuals: reject the coupled response as a whole, across all k points.
1816 50 : IF (.NOT. use_occupation_response .OR. qs_ot_env(1)%gnorm > 0.0_dp) EXIT
1817 12 : use_occupation_response = .FALSE.
1818 : END DO
1819 50 : IF (qs_ot_env(1)%gnorm < 0.0_dp) THEN
1820 0 : logger => cp_get_default_logger()
1821 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1822 : END IF
1823 :
1824 : CALL ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp, &
1825 50 : residual_product=qs_ot_env(1)%gnorm)
1826 :
1827 : ! Detaching a preconditioner leaves the direction workspace in use.
1828 50 : IF (qs_ot_env(1)%use_dx) THEN
1829 38 : IF (do_ks) THEN
1830 88 : DO ispin = 1, nspin
1831 50 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
1832 50 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1833 24 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx_im, -1.0_dp)
1834 : END IF
1835 88 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1836 24 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
1837 24 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
1838 24 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1839 24 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, qs_ot_env(ispin)%rot_mat_gx_im)
1840 24 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx_im, -1.0_dp)
1841 : END IF
1842 : END IF
1843 : END DO
1844 : END IF
1845 38 : IF (do_ener) THEN
1846 36 : DO ispin = 1, nspin
1847 36 : IF (use_occupation_response) THEN
1848 72 : qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_preconditioned_gx
1849 : ELSE
1850 0 : qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_gx
1851 : END IF
1852 : END DO
1853 : END IF
1854 12 : ELSE IF (use_occupation_response) THEN
1855 0 : IF (do_ks) THEN
1856 0 : DO ispin = 1, nspin
1857 0 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_preconditioned_gx)
1858 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1859 0 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_preconditioned_gx_im)
1860 : END IF
1861 : END DO
1862 : END IF
1863 0 : IF (do_ener) THEN
1864 0 : DO ispin = 1, nspin
1865 0 : qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_preconditioned_gx
1866 : END DO
1867 : END IF
1868 : END IF
1869 :
1870 : ! Handling the case of no free variables to optimize
1871 50 : IF (nvariables > 0.0_dp) THEN
1872 48 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
1873 : ELSE
1874 2 : qs_ot_env(1)%gradient = 0.0_dp
1875 : END IF
1876 :
1877 50 : CALL timestop(handle)
1878 :
1879 50 : END SUBROUTINE ot_new_sd_direction
1880 :
1881 : !
1882 : ! creates a new CG direction. Implements Polak-Ribierre variant
1883 : ! using the preconditioner if associated
1884 : ! also updates the gradient for line search
1885 : !
1886 : ! **************************************************************************************************
1887 : !> \brief ...
1888 : !> \param qs_ot_env ...
1889 : !> \param para_env_inter_kp communicator between distributed k-point groups
1890 : ! **************************************************************************************************
1891 23132 : SUBROUTINE ot_new_cg_direction(qs_ot_env, para_env_inter_kp)
1892 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1893 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
1894 :
1895 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_new_cg_direction'
1896 :
1897 : INTEGER :: handle, ispin, nspin
1898 : LOGICAL :: do_ener, do_ks, &
1899 : preceding_response_candidate, &
1900 : use_response_candidate
1901 : REAL(KIND=dp) :: baseline_delta, beta_pr, gnorm_cross, &
1902 : nvariables, test_down, tmp
1903 : TYPE(cp_logger_type), POINTER :: logger
1904 108788 : TYPE(dbcsr_p_type) :: matrix_out(2, SIZE(qs_ot_env))
1905 :
1906 : ! Only the physical low-rank subspace solve can become active. The frozen-H response remains a
1907 : ! shadow when that total projected Hessian is indefinite or unresolved.
1908 :
1909 23132 : CALL timeset(routineN, handle)
1910 :
1911 23132 : nspin = SIZE(qs_ot_env)
1912 23132 : logger => cp_get_default_logger()
1913 23132 : preceding_response_candidate = qs_ot_env(1)%response_candidate_pending
1914 :
1915 23132 : do_ks = qs_ot_env(1)%settings%ks
1916 23132 : do_ener = qs_ot_env(1)%settings%do_ener
1917 23132 : gnorm_cross = 0.0_dp
1918 23132 : IF (do_ks) THEN
1919 51684 : DO ispin = 1, nspin
1920 28552 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
1921 28552 : gnorm_cross = gnorm_cross + tmp
1922 51684 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1923 981 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_gx_old_im, tmp)
1924 981 : gnorm_cross = gnorm_cross + tmp
1925 : END IF
1926 : END DO
1927 23132 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1928 3087 : DO ispin = 1, nspin
1929 1633 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
1930 : ! added 0.5, because we have (antisymmetry) only half the number of variables
1931 1633 : gnorm_cross = gnorm_cross + 0.5_dp*tmp
1932 3087 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1933 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
1934 395 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
1935 395 : gnorm_cross = gnorm_cross + 0.5_dp*tmp
1936 : END IF
1937 : END DO
1938 : END IF
1939 : END IF
1940 23132 : IF (do_ener) THEN
1941 606 : DO ispin = 1, nspin
1942 2578 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
1943 606 : gnorm_cross = gnorm_cross + tmp
1944 : END DO
1945 : END IF
1946 23132 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, gnorm_cross)
1947 :
1948 23132 : qs_ot_env(1)%gnorm = 0.0_dp
1949 23132 : IF (do_ks) THEN
1950 51684 : DO ispin = 1, SIZE(qs_ot_env)
1951 28552 : matrix_out(1, ispin)%matrix => qs_ot_env(ispin)%matrix_gx_old
1952 51684 : matrix_out(2, ispin)%matrix => qs_ot_env(ispin)%matrix_gx_old_im
1953 : END DO
1954 23132 : CALL ot_precondition_gradients(qs_ot_env, matrix_out)
1955 51684 : DO ispin = 1, nspin
1956 28552 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
1957 28552 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1958 51684 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1959 981 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_gx_old_im, tmp)
1960 981 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1961 : END IF
1962 : END DO
1963 23132 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1964 3087 : DO ispin = 1, nspin
1965 1633 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
1966 1633 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
1967 : ! The antisymmetric rotation stores each variable twice.
1968 1633 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
1969 3087 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
1970 395 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old_im, qs_ot_env(ispin)%rot_mat_gx_im)
1971 395 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
1972 395 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
1973 : END IF
1974 : END DO
1975 : END IF
1976 : END IF
1977 23132 : IF (do_ener) THEN
1978 606 : DO ispin = 1, nspin
1979 350 : IF (qs_ot_env(ispin)%settings%occupation_preconditioner) THEN
1980 2578 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_preconditioned_gx
1981 : ELSE
1982 0 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
1983 : END IF
1984 2578 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
1985 606 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1986 : END DO
1987 : END IF
1988 :
1989 23132 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%gnorm)
1990 23132 : IF (.NOT. qs_ot_env(1)%settings%occupation_preconditioner .AND. &
1991 : qs_ot_env(1)%gnorm < 0.0_dp) THEN
1992 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1993 : END IF
1994 :
1995 : ! Occupation-unweighted orbital and energy blocks are useful preconditioner inputs, but only
1996 : ! their coupled direction is relevant for Mermin descent. Fall back as a whole when that
1997 : ! direction is not descending; rejecting individual blocks would destroy the Schur coupling.
1998 23132 : IF (qs_ot_env(1)%settings%occupation_preconditioner .AND. qs_ot_env(1)%gnorm <= 0.0_dp) THEN
1999 2 : qs_ot_env(1)%gnorm = 0.0_dp
2000 2 : IF (do_ks) THEN
2001 2 : CALL ot_precondition_gradients(qs_ot_env, matrix_out, use_occupation_response=.FALSE.)
2002 6 : DO ispin = 1, nspin
2003 4 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
2004 4 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
2005 4 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2006 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
2007 4 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2008 4 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
2009 : END IF
2010 10 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2011 4 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
2012 4 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2013 4 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
2014 4 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2015 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2016 4 : qs_ot_env(ispin)%rot_mat_gx_im)
2017 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2018 4 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2019 4 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
2020 : END IF
2021 : END IF
2022 : END DO
2023 : END IF
2024 2 : IF (do_ener) THEN
2025 6 : DO ispin = 1, nspin
2026 12 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
2027 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + &
2028 : DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
2029 14 : qs_ot_env(ispin)%ener_gx_old)
2030 : END DO
2031 : END IF
2032 2 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%gnorm)
2033 : END IF
2034 : ! Keep convergence reporting tied to the conventional preconditioned residual. The finite
2035 : ! response is an accepted-step candidate, not a redefinition of the physical stopping test.
2036 : CALL ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp, &
2037 23132 : residual_product=qs_ot_env(1)%gnorm)
2038 23132 : baseline_delta = qs_ot_env(1)%delta
2039 : ! Preserve the physical gradients until the stopping test, then retain CG's
2040 : ! preconditioned gradients for its direction and line-search products.
2041 23132 : IF (do_ks .AND. .NOT. qs_ot_env(1)%settings%occupation_preconditioner) THEN
2042 51078 : DO ispin = 1, nspin
2043 51078 : IF (ASSOCIATED(qs_ot_env(ispin)%preconditioner)) THEN
2044 22616 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
2045 22616 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2046 627 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_gx_old_im)
2047 : END IF
2048 : END IF
2049 : END DO
2050 : END IF
2051 23132 : IF (nvariables > 0.0_dp) THEN
2052 23000 : beta_pr = (qs_ot_env(1)%gnorm - gnorm_cross)/qs_ot_env(1)%gnorm_old
2053 : ELSE
2054 132 : beta_pr = 0.0_dp
2055 : END IF
2056 23132 : IF (preceding_response_candidate) beta_pr = 0.0_dp
2057 23132 : IF (cg_history_restart_required( &
2058 : qs_ot_env(1)%settings%occupation_preconditioner, qs_ot_env(1)%etotal, &
2059 2 : qs_ot_env(1)%response_reference_energy)) beta_pr = 0.0_dp
2060 23132 : beta_pr = MAX(beta_pr, 0.0_dp) ! reset to SD
2061 :
2062 23132 : test_down = 0.0_dp
2063 23132 : IF (do_ks) THEN
2064 51684 : DO ispin = 1, nspin
2065 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx_old, &
2066 28552 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2067 28552 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
2068 28552 : test_down = test_down + tmp
2069 28552 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2070 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, qs_ot_env(ispin)%matrix_gx_old_im, &
2071 981 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2072 981 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_dx_im, tmp)
2073 981 : test_down = test_down + tmp
2074 : END IF
2075 80236 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2076 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx_old, &
2077 1633 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2078 1633 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
2079 1633 : test_down = test_down + 0.5_dp*tmp
2080 1633 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2081 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
2082 : qs_ot_env(ispin)%rot_mat_gx_old_im, &
2083 395 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2084 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2085 395 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
2086 395 : test_down = test_down + 0.5_dp*tmp
2087 : END IF
2088 : END IF
2089 : END DO
2090 : END IF
2091 23132 : IF (do_ener) THEN
2092 606 : DO ispin = 1, nspin
2093 : qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - &
2094 2578 : qs_ot_env(ispin)%ener_gx_old
2095 2578 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
2096 606 : test_down = test_down + tmp
2097 : END DO
2098 : END IF
2099 23132 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, test_down)
2100 :
2101 23132 : IF (test_down >= 0.0_dp) THEN ! reset to SD
2102 672 : beta_pr = 0.0_dp
2103 672 : IF (do_ks) THEN
2104 1525 : DO ispin = 1, nspin
2105 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx_old, &
2106 853 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2107 853 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2108 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, qs_ot_env(ispin)%matrix_gx_old_im, &
2109 13 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2110 : END IF
2111 1525 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2112 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
2113 : qs_ot_env(ispin)%rot_mat_gx_old, &
2114 19 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2115 19 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2116 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
2117 : qs_ot_env(ispin)%rot_mat_gx_old_im, &
2118 9 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
2119 : END IF
2120 : END IF
2121 : END DO
2122 : END IF
2123 672 : IF (do_ener) THEN
2124 10 : DO ispin = 1, nspin
2125 : qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - &
2126 58 : qs_ot_env(ispin)%ener_gx_old
2127 : END DO
2128 : END IF
2129 : END IF
2130 :
2131 : CALL ot_try_mermin_response_direction( &
2132 23132 : qs_ot_env, para_env_inter_kp, baseline_delta, test_down, use_response_candidate)
2133 : ! since we change the direction we have to adjust the gradient
2134 23132 : IF (use_response_candidate) THEN
2135 12 : qs_ot_env(1)%OT_METHOD_FULL = "OT CG-R"
2136 12 : qs_ot_env(1)%gradient = test_down
2137 : ELSE
2138 23120 : qs_ot_env(1)%gradient = beta_pr*qs_ot_env(1)%gradient - qs_ot_env(1)%gnorm
2139 : END IF
2140 23132 : qs_ot_env(1)%gnorm_old = qs_ot_env(1)%gnorm
2141 :
2142 23132 : CALL timestop(handle)
2143 :
2144 23132 : END SUBROUTINE ot_new_cg_direction
2145 :
2146 : ! **************************************************************************************************
2147 : !> \brief Returns the smallest shift y <- y + shift*s that meets relative L-BFGS curvature.
2148 : !> \param sy scalar product s.y
2149 : !> \param ss scalar product s.s
2150 : !> \param yy scalar product y.y
2151 : !> \param curvature_tol requested lower bound for (s.y)/sqrt((s.s)(y.y))
2152 : !> \return non-negative finite damping shift, or zero when no valid shift is needed or available
2153 : ! **************************************************************************************************
2154 342 : PURE FUNCTION lbfgs_curvature_damping_shift(sy, ss, yy, curvature_tol) RESULT(damping_shift)
2155 : REAL(KIND=dp), INTENT(IN) :: sy, ss, yy, curvature_tol
2156 : REAL(KIND=dp) :: damping_shift
2157 :
2158 : REAL(KIND=dp) :: curvature_discriminant, target_sy
2159 :
2160 342 : damping_shift = 0.0_dp
2161 : IF (.NOT. ieee_is_finite(sy) .OR. .NOT. ieee_is_finite(ss) .OR. &
2162 : .NOT. ieee_is_finite(yy) .OR. ss <= 0.0_dp .OR. yy <= 0.0_dp .OR. &
2163 342 : curvature_tol < 0.0_dp .OR. curvature_tol >= 1.0_dp) RETURN
2164 :
2165 342 : curvature_discriminant = MAX(0.0_dp, ss*yy - sy*sy)
2166 : target_sy = curvature_tol/SQRT(MAX(TINY(1.0_dp), 1.0_dp - curvature_tol**2))* &
2167 342 : SQRT(curvature_discriminant)
2168 342 : IF (sy <= target_sy) THEN
2169 4 : target_sy = target_sy*(1.0_dp + SQRT(EPSILON(1.0_dp)))
2170 4 : damping_shift = MAX(0.0_dp, (target_sy - sy)/ss)
2171 4 : IF (.NOT. ieee_is_finite(damping_shift)) damping_shift = 0.0_dp
2172 : END IF
2173 : END FUNCTION lbfgs_curvature_damping_shift
2174 :
2175 : ! **************************************************************************************************
2176 : !> \brief Decides whether an L-BFGS history must be discarded after excessive gradient growth.
2177 : !> \param current_gradient_norm_sq squared norm at the current accepted point
2178 : !> \param previous_gradient_norm_sq squared norm at the preceding accepted point
2179 : !> \return true when both norms are finite and the gradient grew by more than a factor of ten
2180 : ! **************************************************************************************************
2181 346 : PURE FUNCTION lbfgs_history_restart_required(current_gradient_norm_sq, &
2182 : previous_gradient_norm_sq) RESULT(restart_history)
2183 : REAL(KIND=dp), INTENT(IN) :: current_gradient_norm_sq, &
2184 : previous_gradient_norm_sq
2185 : LOGICAL :: restart_history
2186 :
2187 : REAL(KIND=dp), PARAMETER :: gradient_restart_factor = 10.0_dp
2188 :
2189 : restart_history = ieee_is_finite(current_gradient_norm_sq) .AND. &
2190 : ieee_is_finite(previous_gradient_norm_sq) .AND. &
2191 : previous_gradient_norm_sq > TINY(1.0_dp) .AND. &
2192 : current_gradient_norm_sq > gradient_restart_factor**2* &
2193 346 : previous_gradient_norm_sq
2194 346 : END FUNCTION lbfgs_history_restart_required
2195 :
2196 : ! **************************************************************************************************
2197 : !> \brief Decides whether L-BFGS must recover from a collapsed accepted line-search step.
2198 : !> \param accepted_step step selected by the preceding line search
2199 : !> \param reference_step configured initial line-search step
2200 : !> \return true for a finite step smaller than one millionth of a valid reference step
2201 : ! **************************************************************************************************
2202 346 : PURE FUNCTION lbfgs_step_restart_required(accepted_step, reference_step) RESULT(restart_step)
2203 : REAL(KIND=dp), INTENT(IN) :: accepted_step, reference_step
2204 : LOGICAL :: restart_step
2205 :
2206 : REAL(KIND=dp), PARAMETER :: step_restart_factor = 1.0E-6_dp
2207 :
2208 : restart_step = ieee_is_finite(accepted_step) .AND. &
2209 : ieee_is_finite(reference_step) .AND. &
2210 : reference_step > 0.0_dp .AND. &
2211 346 : ABS(accepted_step) < step_restart_factor*reference_step
2212 346 : END FUNCTION lbfgs_step_restart_required
2213 :
2214 : ! **************************************************************************************************
2215 : !> \brief Normalize and damp an occupation-response secant against the conventional H0 direction.
2216 : !>
2217 : !> The response direction is first normalized to the product norm of H0*g. If its directional
2218 : !> derivative is smaller than the Powell bound, it is mixed with H0*g until
2219 : !> g^T*p >= 0.2*g^T*H0*g. Equal-norm normalization and convex mixing bound the target direction,
2220 : !> while the curvature bound keeps the inverse-BFGS completion positive and well conditioned.
2221 : !> \param g_dot_response scalar product of g and the raw response direction
2222 : !> \param response_norm_sq squared product norm of the raw response direction
2223 : !> \param g_dot_h0_g scalar product of g and H0*g
2224 : !> \param h0_g_norm_sq squared product norm of H0*g
2225 : !> \param response_scale scale applied to the raw response direction
2226 : !> \param response_weight convex weight of the normalized response direction
2227 : !> \param valid true if finite positive input permits a response secant
2228 : ! **************************************************************************************************
2229 146 : PURE SUBROUTINE lbfgs_response_secant_parameters( &
2230 : g_dot_response, response_norm_sq, g_dot_h0_g, h0_g_norm_sq, &
2231 : response_scale, response_weight, valid)
2232 : REAL(KIND=dp), INTENT(IN) :: g_dot_response, response_norm_sq, &
2233 : g_dot_h0_g, h0_g_norm_sq
2234 : REAL(KIND=dp), INTENT(OUT) :: response_scale, response_weight
2235 : LOGICAL, INTENT(OUT) :: valid
2236 :
2237 : REAL(KIND=dp), PARAMETER :: powell_fraction = 0.2_dp
2238 :
2239 : REAL(KIND=dp) :: scaled_curvature
2240 :
2241 146 : response_scale = 1.0_dp
2242 146 : response_weight = 0.0_dp
2243 : valid = ieee_is_finite(g_dot_response) .AND. &
2244 : ieee_is_finite(response_norm_sq) .AND. &
2245 : ieee_is_finite(g_dot_h0_g) .AND. &
2246 : ieee_is_finite(h0_g_norm_sq) .AND. &
2247 : g_dot_response > 0.0_dp .AND. &
2248 : response_norm_sq > TINY(1.0_dp) .AND. &
2249 : g_dot_h0_g > 0.0_dp .AND. &
2250 146 : h0_g_norm_sq > TINY(1.0_dp)
2251 146 : IF (.NOT. valid) RETURN
2252 :
2253 144 : response_scale = SQRT(h0_g_norm_sq/response_norm_sq)
2254 144 : scaled_curvature = response_scale*g_dot_response
2255 : IF (.NOT. ieee_is_finite(response_scale) .OR. &
2256 144 : .NOT. ieee_is_finite(scaled_curvature) .OR. scaled_curvature <= 0.0_dp) THEN
2257 0 : response_scale = 1.0_dp
2258 0 : valid = .FALSE.
2259 0 : RETURN
2260 : END IF
2261 :
2262 144 : response_weight = 1.0_dp
2263 144 : IF (scaled_curvature < powell_fraction*g_dot_h0_g) THEN
2264 : response_weight = (1.0_dp - powell_fraction)*g_dot_h0_g/ &
2265 66 : (g_dot_h0_g - scaled_curvature)
2266 : END IF
2267 : END SUBROUTINE lbfgs_response_secant_parameters
2268 :
2269 : ! **************************************************************************************************
2270 : !> \brief Apply the L-BFGS initial inverse Hessian in the OT product space.
2271 : !>
2272 : !> The conventional OT preconditioner defines the base operator H0. If occupation response is
2273 : !> enabled, a norm-bounded Powell-damped response direction p defines a positive inverse-BFGS
2274 : !> secant completion. It maps the current physical gradient g to p while remaining a linear
2275 : !> operator for the arbitrary vector q produced by the first L-BFGS loop:
2276 : !>
2277 : !> H = (I - rho*p*g^T)*H0*(I - rho*g*p^T) + rho*p*p^T,
2278 : !>
2279 : !> where rho = 1/(g^T*p). The m+1 history slot is scratch space for p; callers restore the current
2280 : !> physical gradient there after this routine returns.
2281 : !> \param qs_ot_env OT environments for all local spin/k-point channels
2282 : !> \param gamma scalar H0 scale for orbital and energy variables without an OT preconditioner
2283 : !> \param gamma_rotation H0 scale for rotation variables
2284 : !> \param scratch_index history slot used to store p
2285 : !> \param para_env_inter_kp communicator between distributed k-point groups
2286 : ! **************************************************************************************************
2287 600 : SUBROUTINE ot_apply_lbfgs_initial_inverse( &
2288 : qs_ot_env, gamma, gamma_rotation, scratch_index, para_env_inter_kp)
2289 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
2290 : REAL(KIND=dp), INTENT(IN) :: gamma, gamma_rotation
2291 : INTEGER, INTENT(IN) :: scratch_index
2292 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
2293 :
2294 : INTEGER :: ispin
2295 : LOGICAL :: valid_response
2296 : REAL(KIND=dp) :: coeff_p, coeff_v, g_dot_h0_g, g_dot_h0_q, g_dot_p, h0_g_norm_sq, p_dot_q, &
2297 : p_norm_sq, response_scale, response_weight, rho, tmp
2298 : TYPE(dbcsr_type), POINTER :: scratch_im
2299 :
2300 : ! Apply the base operator to q, which is stored in the *_gx_old work arrays.
2301 1154 : DO ispin = 1, SIZE(qs_ot_env)
2302 : CALL ot_apply_channel_preconditioner(qs_ot_env(ispin), &
2303 : qs_ot_env(ispin)%matrix_gx_old, &
2304 : qs_ot_env(ispin)%matrix_dx, &
2305 : qs_ot_env(ispin)%matrix_gx_old_im, &
2306 694 : qs_ot_env(ispin)%matrix_dx_im, identity_scale=gamma)
2307 694 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2308 566 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx_old)
2309 566 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, gamma_rotation)
2310 566 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2311 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, &
2312 566 : qs_ot_env(ispin)%rot_mat_gx_old_im)
2313 566 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx_im, gamma_rotation)
2314 : END IF
2315 : END IF
2316 1154 : IF (qs_ot_env(1)%settings%do_ener) THEN
2317 2612 : qs_ot_env(ispin)%ener_dx = gamma*qs_ot_env(ispin)%ener_gx_old
2318 : END IF
2319 : END DO
2320 460 : IF (.NOT. qs_ot_env(1)%settings%occupation_preconditioner) RETURN
2321 :
2322 : ! Build the desired occupation-preconditioned product direction p in the scratch slot.
2323 344 : DO ispin = 1, SIZE(qs_ot_env)
2324 204 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%matrix_preconditioned_gx))
2325 204 : NULLIFY (scratch_im)
2326 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2327 204 : scratch_im => qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix
2328 : END IF
2329 : CALL ot_apply_channel_preconditioner(qs_ot_env(ispin), &
2330 : qs_ot_env(ispin)%matrix_preconditioned_gx, &
2331 : qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2332 : qs_ot_env(ispin)%matrix_preconditioned_gx_im, &
2333 204 : scratch_im)
2334 204 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2335 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2336 204 : qs_ot_env(ispin)%rot_mat_gx)
2337 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2338 204 : gamma_rotation)
2339 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2340 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2341 204 : qs_ot_env(ispin)%rot_mat_gx_im)
2342 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2343 204 : gamma_rotation)
2344 : END IF
2345 : END IF
2346 344 : IF (qs_ot_env(1)%settings%do_ener) THEN
2347 204 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%ener_preconditioned_gx))
2348 : qs_ot_env(ispin)%ener_h_e(scratch_index, :) = &
2349 1140 : qs_ot_env(ispin)%ener_preconditioned_gx
2350 : END IF
2351 : END DO
2352 :
2353 140 : g_dot_h0_q = 0.0_dp
2354 140 : g_dot_p = 0.0_dp
2355 140 : p_dot_q = 0.0_dp
2356 140 : p_norm_sq = 0.0_dp
2357 344 : DO ispin = 1, SIZE(qs_ot_env)
2358 204 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
2359 204 : g_dot_h0_q = g_dot_h0_q + tmp
2360 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
2361 204 : qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, tmp)
2362 204 : g_dot_p = g_dot_p + tmp
2363 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2364 204 : qs_ot_env(ispin)%matrix_gx_old, tmp)
2365 204 : p_dot_q = p_dot_q + tmp
2366 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2367 204 : qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, tmp)
2368 204 : p_norm_sq = p_norm_sq + tmp
2369 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2370 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
2371 204 : qs_ot_env(ispin)%matrix_dx_im, tmp)
2372 204 : g_dot_h0_q = g_dot_h0_q + tmp
2373 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
2374 204 : qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, tmp)
2375 204 : g_dot_p = g_dot_p + tmp
2376 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, &
2377 204 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2378 204 : p_dot_q = p_dot_q + tmp
2379 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, &
2380 204 : qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, tmp)
2381 204 : p_norm_sq = p_norm_sq + tmp
2382 : END IF
2383 204 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2384 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
2385 204 : qs_ot_env(ispin)%rot_mat_dx, tmp)
2386 204 : g_dot_h0_q = g_dot_h0_q + 0.5_dp*tmp
2387 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
2388 204 : qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, tmp)
2389 204 : g_dot_p = g_dot_p + 0.5_dp*tmp
2390 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2391 204 : qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2392 204 : p_dot_q = p_dot_q + 0.5_dp*tmp
2393 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2394 204 : qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, tmp)
2395 204 : p_norm_sq = p_norm_sq + 0.5_dp*tmp
2396 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2397 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2398 204 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
2399 204 : g_dot_h0_q = g_dot_h0_q + 0.5_dp*tmp
2400 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2401 204 : qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, tmp)
2402 204 : g_dot_p = g_dot_p + 0.5_dp*tmp
2403 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2404 204 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2405 204 : p_dot_q = p_dot_q + 0.5_dp*tmp
2406 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2407 204 : qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, tmp)
2408 204 : p_norm_sq = p_norm_sq + 0.5_dp*tmp
2409 : END IF
2410 : END IF
2411 548 : IF (qs_ot_env(1)%settings%do_ener) THEN
2412 : g_dot_h0_q = g_dot_h0_q + &
2413 1140 : DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
2414 : g_dot_p = g_dot_p + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
2415 1140 : qs_ot_env(ispin)%ener_h_e(scratch_index, :))
2416 : p_dot_q = p_dot_q + DOT_PRODUCT(qs_ot_env(ispin)%ener_h_e(scratch_index, :), &
2417 1140 : qs_ot_env(ispin)%ener_gx_old)
2418 : p_norm_sq = p_norm_sq + DOT_PRODUCT( &
2419 : qs_ot_env(ispin)%ener_h_e(scratch_index, :), &
2420 1140 : qs_ot_env(ispin)%ener_h_e(scratch_index, :))
2421 : END IF
2422 : END DO
2423 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_h0_q)
2424 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_p)
2425 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, p_dot_q)
2426 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, p_norm_sq)
2427 :
2428 : ! Overwrite q with H0*g; q is no longer needed after p^T*q has been formed.
2429 140 : g_dot_h0_g = 0.0_dp
2430 140 : h0_g_norm_sq = 0.0_dp
2431 344 : DO ispin = 1, SIZE(qs_ot_env)
2432 : CALL ot_apply_channel_preconditioner(qs_ot_env(ispin), &
2433 : qs_ot_env(ispin)%matrix_gx, &
2434 : qs_ot_env(ispin)%matrix_gx_old, &
2435 : qs_ot_env(ispin)%matrix_gx_im, &
2436 204 : qs_ot_env(ispin)%matrix_gx_old_im, identity_scale=gamma)
2437 204 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
2438 204 : g_dot_h0_g = g_dot_h0_g + tmp
2439 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_old, &
2440 204 : qs_ot_env(ispin)%matrix_gx_old, tmp)
2441 204 : h0_g_norm_sq = h0_g_norm_sq + tmp
2442 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2443 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
2444 204 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2445 204 : g_dot_h0_g = g_dot_h0_g + tmp
2446 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_old_im, &
2447 204 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2448 204 : h0_g_norm_sq = h0_g_norm_sq + tmp
2449 : END IF
2450 204 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2451 204 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
2452 204 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_gx_old, gamma_rotation)
2453 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
2454 204 : qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2455 204 : g_dot_h0_g = g_dot_h0_g + 0.5_dp*tmp
2456 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_old, &
2457 204 : qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2458 204 : h0_g_norm_sq = h0_g_norm_sq + 0.5_dp*tmp
2459 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2460 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2461 204 : qs_ot_env(ispin)%rot_mat_gx_im)
2462 204 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_gx_old_im, gamma_rotation)
2463 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2464 204 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2465 204 : g_dot_h0_g = g_dot_h0_g + 0.5_dp*tmp
2466 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2467 204 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2468 204 : h0_g_norm_sq = h0_g_norm_sq + 0.5_dp*tmp
2469 : END IF
2470 : END IF
2471 548 : IF (qs_ot_env(1)%settings%do_ener) THEN
2472 1140 : qs_ot_env(ispin)%ener_gx_old = gamma*qs_ot_env(ispin)%ener_gx
2473 : g_dot_h0_g = g_dot_h0_g + &
2474 : DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
2475 1140 : qs_ot_env(ispin)%ener_gx_old)
2476 : h0_g_norm_sq = h0_g_norm_sq + &
2477 : DOT_PRODUCT(qs_ot_env(ispin)%ener_gx_old, &
2478 1140 : qs_ot_env(ispin)%ener_gx_old)
2479 : END IF
2480 : END DO
2481 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_h0_g)
2482 140 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, h0_g_norm_sq)
2483 :
2484 : CALL lbfgs_response_secant_parameters( &
2485 : g_dot_p, p_norm_sq, g_dot_h0_g, h0_g_norm_sq, &
2486 140 : response_scale, response_weight, valid_response)
2487 140 : IF (.NOT. valid_response) RETURN
2488 :
2489 : ! Replace p by the bounded, damped response target and update its scalar products.
2490 344 : DO ispin = 1, SIZE(qs_ot_env)
2491 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2492 204 : response_weight*response_scale)
2493 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2494 : qs_ot_env(ispin)%matrix_gx_old, &
2495 204 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp - response_weight)
2496 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2497 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, &
2498 204 : response_weight*response_scale)
2499 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, &
2500 : qs_ot_env(ispin)%matrix_gx_old_im, &
2501 204 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp - response_weight)
2502 : END IF
2503 204 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2504 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2505 204 : response_weight*response_scale)
2506 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2507 : qs_ot_env(ispin)%rot_mat_gx_old, &
2508 204 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp - response_weight)
2509 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2510 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2511 204 : response_weight*response_scale)
2512 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2513 : qs_ot_env(ispin)%rot_mat_gx_old_im, &
2514 204 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp - response_weight)
2515 : END IF
2516 : END IF
2517 344 : IF (qs_ot_env(1)%settings%do_ener) THEN
2518 : qs_ot_env(ispin)%ener_h_e(scratch_index, :) = &
2519 : response_weight*response_scale*qs_ot_env(ispin)%ener_h_e(scratch_index, :) + &
2520 1140 : (1.0_dp - response_weight)*qs_ot_env(ispin)%ener_gx_old
2521 : END IF
2522 : END DO
2523 : g_dot_p = response_weight*response_scale*g_dot_p + &
2524 140 : (1.0_dp - response_weight)*g_dot_h0_g
2525 : p_dot_q = response_weight*response_scale*p_dot_q + &
2526 140 : (1.0_dp - response_weight)*g_dot_h0_q
2527 :
2528 140 : rho = 1.0_dp/g_dot_p
2529 140 : coeff_p = -rho*g_dot_h0_q + p_dot_q*(rho + rho*rho*g_dot_h0_g)
2530 140 : coeff_v = -rho*p_dot_q
2531 344 : DO ispin = 1, SIZE(qs_ot_env)
2532 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, &
2533 : qs_ot_env(ispin)%matrix_h_e(scratch_index)%matrix, &
2534 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_p)
2535 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx_old, &
2536 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_v)
2537 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2538 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, &
2539 : qs_ot_env(ispin)%matrix_h_e_im(scratch_index)%matrix, &
2540 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_p)
2541 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, &
2542 : qs_ot_env(ispin)%matrix_gx_old_im, &
2543 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_v)
2544 : END IF
2545 204 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2546 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
2547 : qs_ot_env(ispin)%rot_mat_h_e(scratch_index)%matrix, &
2548 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_p)
2549 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx_old, &
2550 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_v)
2551 204 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2552 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
2553 : qs_ot_env(ispin)%rot_mat_h_e_im(scratch_index)%matrix, &
2554 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_p)
2555 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
2556 : qs_ot_env(ispin)%rot_mat_gx_old_im, &
2557 204 : alpha_scalar=1.0_dp, beta_scalar=coeff_v)
2558 : END IF
2559 : END IF
2560 344 : IF (qs_ot_env(1)%settings%do_ener) THEN
2561 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_dx + &
2562 : coeff_p*qs_ot_env(ispin)%ener_h_e(scratch_index, :) + &
2563 1140 : coeff_v*qs_ot_env(ispin)%ener_gx_old
2564 : END IF
2565 : END DO
2566 : END SUBROUTINE ot_apply_lbfgs_initial_inverse
2567 :
2568 : ! **************************************************************************************************
2569 : !> \brief Builds an L-BFGS direction in the fixed OT product chart.
2570 : !>
2571 : !> The selected OT preconditioner is used as the initial inverse Hessian. History slot m+1 stores
2572 : !> the previous accepted point and gradient; slots 1:m store the circular (s,y) secant history.
2573 : !> Orbital rotations are included with the antisymmetric-matrix metric used by the other OT
2574 : !> minimizers. Weak curvature can be regularized along the step, while rejected updates leave
2575 : !> older valid history intact. Non-descent directions and collapsed line-search steps reset the
2576 : !> history; the latter also restore the configured initial step.
2577 : !> \param qs_ot_env OT environments for all local spin/k-point channels
2578 : !> \param para_env_inter_kp communicator between distributed k-point groups
2579 : ! **************************************************************************************************
2580 460 : SUBROUTINE ot_new_lbfgs_direction(qs_ot_env, para_env_inter_kp)
2581 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
2582 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
2583 :
2584 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_new_lbfgs_direction'
2585 : REAL(KIND=dp), PARAMETER :: rotation_scale_max = 1.0E3_dp, &
2586 : rotation_scale_min = 1.0E-3_dp
2587 :
2588 : INTEGER :: handle, i, ispin, j, m_history, newest, &
2589 : nhistory, nspin
2590 460 : INTEGER, ALLOCATABLE, DIMENSION(:) :: history_index
2591 : LOGICAL :: do_ener, preceding_response_candidate, &
2592 : restart_history, restart_step, &
2593 : use_response_candidate
2594 : REAL(KIND=dp) :: beta, current_gradient_norm_sq, curvature_tol, damping_shift, g_dot_z, &
2595 : gamma, gamma_rotation, nvariables, previous_gradient_norm_sq, ss, ss_rotation, stq, sy, &
2596 : sy_rotation, test_down, tmp, yy, yy_newest, yy_rotation
2597 460 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, rho
2598 :
2599 460 : CALL timeset(routineN, handle)
2600 :
2601 460 : IF (.NOT. qs_ot_env(1)%settings%ks) THEN
2602 0 : CPABORT("MINIMIZER LBFGS currently requires OT orbital variables")
2603 : END IF
2604 460 : nspin = SIZE(qs_ot_env)
2605 460 : do_ener = qs_ot_env(1)%settings%do_ener
2606 460 : m_history = qs_ot_env(1)%settings%diis_m
2607 460 : curvature_tol = qs_ot_env(1)%settings%lbfgs_curvature_tol
2608 460 : restart_history = .FALSE.
2609 460 : restart_step = .FALSE.
2610 460 : preceding_response_candidate = qs_ot_env(1)%response_candidate_pending
2611 :
2612 : ! A local line search can accept a nearly isoenergetic point whose raw gradient is much
2613 : ! larger than at the preceding point. CG naturally resets after such a loss of conjugacy,
2614 : ! whereas retaining an L-BFGS model can amplify it. Discard only the secant model and use
2615 : ! the physical initial inverse Hessian for the current direction. Likewise, restore the
2616 : ! configured initial step if repeated interpolation has collapsed the accepted step by more
2617 : ! than six orders of magnitude without satisfying the SCF convergence criterion.
2618 460 : IF (qs_ot_env(1)%diis_iter > 0) THEN
2619 342 : current_gradient_norm_sq = 0.0_dp
2620 342 : previous_gradient_norm_sq = 0.0_dp
2621 886 : DO ispin = 1, nspin
2622 544 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
2623 544 : current_gradient_norm_sq = current_gradient_norm_sq + tmp
2624 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e(m_history + 1)%matrix, &
2625 544 : qs_ot_env(ispin)%matrix_h_e(m_history + 1)%matrix, tmp)
2626 544 : previous_gradient_norm_sq = previous_gradient_norm_sq + tmp
2627 544 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2628 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
2629 544 : qs_ot_env(ispin)%matrix_gx_im, tmp)
2630 544 : current_gradient_norm_sq = current_gradient_norm_sq + tmp
2631 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e_im(m_history + 1)%matrix, &
2632 544 : qs_ot_env(ispin)%matrix_h_e_im(m_history + 1)%matrix, tmp)
2633 544 : previous_gradient_norm_sq = previous_gradient_norm_sq + tmp
2634 : END IF
2635 544 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2636 504 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
2637 504 : current_gradient_norm_sq = current_gradient_norm_sq + 0.5_dp*tmp
2638 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e(m_history + 1)%matrix, &
2639 504 : qs_ot_env(ispin)%rot_mat_h_e(m_history + 1)%matrix, tmp)
2640 504 : previous_gradient_norm_sq = previous_gradient_norm_sq + 0.5_dp*tmp
2641 504 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2642 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
2643 504 : qs_ot_env(ispin)%rot_mat_gx_im, tmp)
2644 504 : current_gradient_norm_sq = current_gradient_norm_sq + 0.5_dp*tmp
2645 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e_im(m_history + 1)%matrix, &
2646 504 : qs_ot_env(ispin)%rot_mat_h_e_im(m_history + 1)%matrix, tmp)
2647 504 : previous_gradient_norm_sq = previous_gradient_norm_sq + 0.5_dp*tmp
2648 : END IF
2649 : END IF
2650 1430 : IF (do_ener) THEN
2651 2330 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
2652 242 : current_gradient_norm_sq = current_gradient_norm_sq + tmp
2653 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_e(m_history + 1, :), &
2654 2330 : qs_ot_env(ispin)%ener_h_e(m_history + 1, :))
2655 242 : previous_gradient_norm_sq = previous_gradient_norm_sq + tmp
2656 : END IF
2657 : END DO
2658 342 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, current_gradient_norm_sq)
2659 342 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, previous_gradient_norm_sq)
2660 : restart_step = lbfgs_step_restart_required(qs_ot_env(1)%ds_min, &
2661 342 : qs_ot_env(1)%settings%ds_min)
2662 : restart_history = lbfgs_history_restart_required(current_gradient_norm_sq, &
2663 : previous_gradient_norm_sq) .OR. &
2664 342 : restart_step .OR. preceding_response_candidate
2665 : IF (restart_history) THEN
2666 6 : qs_ot_env(1)%diis_iter = 1
2667 6 : IF (restart_step) qs_ot_env(1)%ds_min = qs_ot_env(1)%settings%ds_min
2668 6 : qs_ot_env(1)%OT_METHOD_FULL = "OT L-RST"
2669 : END IF
2670 : END IF
2671 :
2672 : ! Form a candidate secant pair in the direction workspaces relative to the preceding
2673 : ! accepted line-search point. Commit it to the circular history only after the curvature
2674 : ! test, so a rejected candidate cannot overwrite the oldest valid pair in a full buffer.
2675 460 : IF (qs_ot_env(1)%diis_iter > 0 .AND. .NOT. restart_history) THEN
2676 336 : j = MOD(qs_ot_env(1)%diis_iter - 1, m_history) + 1
2677 336 : ss = 0.0_dp
2678 336 : sy = 0.0_dp
2679 336 : sy_rotation = 0.0_dp
2680 336 : ss_rotation = 0.0_dp
2681 336 : yy = 0.0_dp
2682 336 : yy_rotation = 0.0_dp
2683 874 : DO ispin = 1, nspin
2684 538 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_x)
2685 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, &
2686 : qs_ot_env(ispin)%matrix_h_x(m_history + 1)%matrix, &
2687 538 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2688 538 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx)
2689 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old, &
2690 : qs_ot_env(ispin)%matrix_h_e(m_history + 1)%matrix, &
2691 538 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2692 538 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_dx, tmp)
2693 538 : ss = ss + tmp
2694 538 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx_old, tmp)
2695 538 : sy = sy + tmp
2696 538 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx_old, tmp)
2697 538 : yy = yy + tmp
2698 538 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2699 538 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx_im, qs_ot_env(ispin)%matrix_x_im)
2700 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, &
2701 : qs_ot_env(ispin)%matrix_h_x_im(m_history + 1)%matrix, &
2702 538 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2703 538 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old_im, qs_ot_env(ispin)%matrix_gx_im)
2704 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old_im, &
2705 : qs_ot_env(ispin)%matrix_h_e_im(m_history + 1)%matrix, &
2706 538 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2707 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_dx_im, &
2708 538 : qs_ot_env(ispin)%matrix_dx_im, tmp)
2709 538 : ss = ss + tmp
2710 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_dx_im, &
2711 538 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2712 538 : sy = sy + tmp
2713 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_old_im, &
2714 538 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2715 538 : yy = yy + tmp
2716 : END IF
2717 538 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2718 498 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_x)
2719 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
2720 : qs_ot_env(ispin)%rot_mat_h_x(m_history + 1)%matrix, &
2721 498 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2722 498 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
2723 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old, &
2724 : qs_ot_env(ispin)%rot_mat_h_e(m_history + 1)%matrix, &
2725 498 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2726 498 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_dx, tmp)
2727 498 : ss = ss + 0.5_dp*tmp
2728 498 : ss_rotation = ss_rotation + 0.5_dp*tmp
2729 498 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2730 498 : sy = sy + 0.5_dp*tmp
2731 498 : sy_rotation = sy_rotation + 0.5_dp*tmp
2732 498 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2733 498 : yy = yy + 0.5_dp*tmp
2734 498 : yy_rotation = yy_rotation + 0.5_dp*tmp
2735 498 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2736 498 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, qs_ot_env(ispin)%rot_mat_x_im)
2737 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
2738 : qs_ot_env(ispin)%rot_mat_h_x_im(m_history + 1)%matrix, &
2739 498 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2740 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2741 498 : qs_ot_env(ispin)%rot_mat_gx_im)
2742 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2743 : qs_ot_env(ispin)%rot_mat_h_e_im(m_history + 1)%matrix, &
2744 498 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2745 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_dx_im, &
2746 498 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
2747 498 : ss = ss + 0.5_dp*tmp
2748 498 : ss_rotation = ss_rotation + 0.5_dp*tmp
2749 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_dx_im, &
2750 498 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2751 498 : sy = sy + 0.5_dp*tmp
2752 498 : sy_rotation = sy_rotation + 0.5_dp*tmp
2753 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2754 498 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2755 498 : yy = yy + 0.5_dp*tmp
2756 498 : yy_rotation = yy_rotation + 0.5_dp*tmp
2757 : END IF
2758 : END IF
2759 1412 : IF (do_ener) THEN
2760 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_x - &
2761 2288 : qs_ot_env(ispin)%ener_h_x(m_history + 1, :)
2762 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx - &
2763 2288 : qs_ot_env(ispin)%ener_h_e(m_history + 1, :)
2764 2288 : ss = ss + DOT_PRODUCT(qs_ot_env(ispin)%ener_dx, qs_ot_env(ispin)%ener_dx)
2765 2288 : sy = sy + DOT_PRODUCT(qs_ot_env(ispin)%ener_dx, qs_ot_env(ispin)%ener_gx_old)
2766 2288 : yy = yy + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx_old, qs_ot_env(ispin)%ener_gx_old)
2767 : END IF
2768 : END DO
2769 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, ss)
2770 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, sy)
2771 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, yy)
2772 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, ss_rotation)
2773 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, sy_rotation)
2774 336 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, yy_rotation)
2775 :
2776 : ! Shift y along s just enough to satisfy the requested relative curvature. With the
2777 : ! product-chart metric this preserves orbital-rotation covariance. The target below is
2778 : ! obtained by solving (s.y)^2 >= tol^2 (s.s)(y.y) for y <- y + shift*s.
2779 336 : damping_shift = 0.0_dp
2780 336 : IF (qs_ot_env(1)%settings%lbfgs_damping) THEN
2781 336 : damping_shift = lbfgs_curvature_damping_shift(sy, ss, yy, curvature_tol)
2782 : END IF
2783 336 : IF (damping_shift > 0.0_dp) THEN
2784 0 : DO ispin = 1, nspin
2785 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old, &
2786 : qs_ot_env(ispin)%matrix_dx, &
2787 0 : alpha_scalar=1.0_dp, beta_scalar=damping_shift)
2788 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2789 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old_im, &
2790 : qs_ot_env(ispin)%matrix_dx_im, &
2791 0 : alpha_scalar=1.0_dp, beta_scalar=damping_shift)
2792 : END IF
2793 0 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2794 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old, &
2795 : qs_ot_env(ispin)%rot_mat_dx, &
2796 0 : alpha_scalar=1.0_dp, beta_scalar=damping_shift)
2797 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2798 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2799 : qs_ot_env(ispin)%rot_mat_dx_im, &
2800 0 : alpha_scalar=1.0_dp, beta_scalar=damping_shift)
2801 : END IF
2802 : END IF
2803 0 : IF (do_ener) THEN
2804 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx_old + &
2805 0 : damping_shift*qs_ot_env(ispin)%ener_dx
2806 : END IF
2807 : END DO
2808 0 : yy = yy + 2.0_dp*damping_shift*sy + damping_shift**2*ss
2809 : yy_rotation = yy_rotation + 2.0_dp*damping_shift*sy_rotation + &
2810 0 : damping_shift**2*ss_rotation
2811 0 : sy = sy + damping_shift*ss
2812 0 : sy_rotation = sy_rotation + damping_shift*ss_rotation
2813 : END IF
2814 : IF (ieee_is_finite(sy) .AND. ieee_is_finite(ss) .AND. ieee_is_finite(yy) .AND. &
2815 336 : ss > 0.0_dp .AND. yy > 0.0_dp .AND. sy > curvature_tol*SQRT(ss*yy)) THEN
2816 874 : DO ispin = 1, nspin
2817 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, &
2818 538 : qs_ot_env(ispin)%matrix_dx)
2819 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
2820 538 : qs_ot_env(ispin)%matrix_gx_old)
2821 538 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2822 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x_im(j)%matrix, &
2823 538 : qs_ot_env(ispin)%matrix_dx_im)
2824 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_e_im(j)%matrix, &
2825 538 : qs_ot_env(ispin)%matrix_gx_old_im)
2826 : END IF
2827 538 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2828 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, &
2829 498 : qs_ot_env(ispin)%rot_mat_dx)
2830 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, &
2831 498 : qs_ot_env(ispin)%rot_mat_gx_old)
2832 498 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2833 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x_im(j)%matrix, &
2834 498 : qs_ot_env(ispin)%rot_mat_dx_im)
2835 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e_im(j)%matrix, &
2836 498 : qs_ot_env(ispin)%rot_mat_gx_old_im)
2837 : END IF
2838 : END IF
2839 874 : IF (do_ener) THEN
2840 2288 : qs_ot_env(ispin)%ener_h_x(j, :) = qs_ot_env(ispin)%ener_dx
2841 2288 : qs_ot_env(ispin)%ener_h_e(j, :) = qs_ot_env(ispin)%ener_gx_old
2842 : END IF
2843 : END DO
2844 336 : qs_ot_env(1)%lbfgs_rho(j) = 1.0_dp/sy
2845 336 : qs_ot_env(1)%lbfgs_yy(j) = yy
2846 336 : qs_ot_env(1)%lbfgs_sy_rotation(j) = sy_rotation
2847 336 : qs_ot_env(1)%lbfgs_yy_rotation(j) = yy_rotation
2848 336 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
2849 : ELSE
2850 : ! Skip this pair. Older pairs remain valid in the fixed OT chart and are retained;
2851 : ! the descent safeguard below still clears all history if their recursion fails.
2852 0 : qs_ot_env(1)%OT_METHOD_FULL = "OT LSKIP"
2853 : END IF
2854 124 : ELSE IF (qs_ot_env(1)%diis_iter <= 0) THEN
2855 118 : qs_ot_env(1)%diis_iter = 1
2856 : END IF
2857 :
2858 1154 : DO ispin = 1, nspin
2859 694 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx)
2860 694 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2861 694 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old_im, qs_ot_env(ispin)%matrix_gx_im)
2862 : END IF
2863 694 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2864 566 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
2865 566 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2866 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2867 566 : qs_ot_env(ispin)%rot_mat_gx_im)
2868 : END IF
2869 : END IF
2870 3498 : IF (do_ener) qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
2871 : END DO
2872 :
2873 460 : nhistory = MIN(qs_ot_env(1)%diis_iter - 1, m_history)
2874 : ALLOCATE (alpha(MAX(1, nhistory)), rho(MAX(1, nhistory)), &
2875 2760 : history_index(MAX(1, nhistory)))
2876 460 : gamma_rotation = 1.0_dp
2877 460 : yy_newest = 0.0_dp
2878 :
2879 : ! First loop, newest to oldest: q <- V^T g.
2880 460 : IF (nhistory > 0) THEN
2881 336 : newest = MOD(qs_ot_env(1)%diis_iter - 2, m_history) + 1
2882 336 : yy_newest = qs_ot_env(1)%lbfgs_yy(newest)
2883 : IF (qs_ot_env(1)%settings%do_rotation .AND. &
2884 : ieee_is_finite(qs_ot_env(1)%lbfgs_sy_rotation(newest)) .AND. &
2885 : ieee_is_finite(qs_ot_env(1)%lbfgs_yy_rotation(newest)) .AND. &
2886 336 : qs_ot_env(1)%lbfgs_sy_rotation(newest) > 0.0_dp .AND. &
2887 : qs_ot_env(1)%lbfgs_yy_rotation(newest) > 0.0_dp) THEN
2888 : gamma_rotation = MAX(rotation_scale_min, MIN(rotation_scale_max, &
2889 : qs_ot_env(1)%lbfgs_sy_rotation(newest)/ &
2890 308 : qs_ot_env(1)%lbfgs_yy_rotation(newest)))
2891 : END IF
2892 1842 : DO i = 1, nhistory
2893 1506 : history_index(i) = MODULO(newest - i, m_history) + 1
2894 1506 : stq = 0.0_dp
2895 3754 : DO ispin = 1, nspin
2896 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_x(history_index(i))%matrix, &
2897 2248 : qs_ot_env(ispin)%matrix_gx_old, tmp)
2898 2248 : stq = stq + tmp
2899 2248 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2900 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_x_im(history_index(i))%matrix, &
2901 2248 : qs_ot_env(ispin)%matrix_gx_old_im, tmp)
2902 2248 : stq = stq + tmp
2903 : END IF
2904 2248 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2905 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_x(history_index(i))%matrix, &
2906 2152 : qs_ot_env(ispin)%rot_mat_gx_old, tmp)
2907 2152 : stq = stq + 0.5_dp*tmp
2908 2152 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2909 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_x_im(history_index(i))%matrix, &
2910 2152 : qs_ot_env(ispin)%rot_mat_gx_old_im, tmp)
2911 2152 : stq = stq + 0.5_dp*tmp
2912 : END IF
2913 : END IF
2914 6002 : IF (do_ener) THEN
2915 : stq = stq + DOT_PRODUCT(qs_ot_env(ispin)%ener_h_x(history_index(i), :), &
2916 10278 : qs_ot_env(ispin)%ener_gx_old)
2917 : END IF
2918 : END DO
2919 1506 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, stq)
2920 1506 : rho(i) = qs_ot_env(1)%lbfgs_rho(history_index(i))
2921 1506 : alpha(i) = rho(i)*stq
2922 4090 : DO ispin = 1, nspin
2923 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old, &
2924 : qs_ot_env(ispin)%matrix_h_e(history_index(i))%matrix, &
2925 2248 : alpha_scalar=1.0_dp, beta_scalar=-alpha(i))
2926 2248 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2927 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_gx_old_im, &
2928 : qs_ot_env(ispin)%matrix_h_e_im(history_index(i))%matrix, &
2929 2248 : alpha_scalar=1.0_dp, beta_scalar=-alpha(i))
2930 : END IF
2931 2248 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2932 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old, &
2933 : qs_ot_env(ispin)%rot_mat_h_e(history_index(i))%matrix, &
2934 2152 : alpha_scalar=1.0_dp, beta_scalar=-alpha(i))
2935 2152 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2936 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_gx_old_im, &
2937 : qs_ot_env(ispin)%rot_mat_h_e_im(history_index(i))%matrix, &
2938 2152 : alpha_scalar=1.0_dp, beta_scalar=-alpha(i))
2939 : END IF
2940 : END IF
2941 3754 : IF (do_ener) THEN
2942 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx_old - &
2943 10278 : alpha(i)*qs_ot_env(ispin)%ener_h_e(history_index(i), :)
2944 : END IF
2945 : END DO
2946 : END DO
2947 : END IF
2948 :
2949 : ! Apply the conventional initial inverse Hessian and, when requested, its positive linear
2950 : ! occupation-response completion to the arbitrary vector from the first loop.
2951 124 : gamma = 1.0_dp
2952 336 : IF (nhistory > 0 .AND. yy_newest > 0.0_dp) gamma = 1.0_dp/(rho(1)*yy_newest)
2953 : CALL ot_apply_lbfgs_initial_inverse( &
2954 460 : qs_ot_env, gamma, gamma_rotation, m_history + 1, para_env_inter_kp)
2955 :
2956 : ! Real-orbital convergence uses the conventional coupled residual below.
2957 460 : IF (.NOT. qs_ot_env(1)%has_complex_kpoint_state) THEN
2958 : CALL ot_build_history_residual(qs_ot_env, m_history + 1, 1.0_dp, &
2959 0 : para_env_inter_kp, qs_ot_env(1)%gnorm)
2960 0 : IF (.NOT. ieee_is_finite(qs_ot_env(1)%gnorm) .OR. qs_ot_env(1)%gnorm <= 0.0_dp) THEN
2961 : ! Match the identity-metric descent fallback without using the history-based direction.
2962 : CALL ot_build_history_residual_once(qs_ot_env, m_history + 1, 1.0_dp, .FALSE., &
2963 0 : qs_ot_env(1)%gnorm, use_preconditioner=.FALSE.)
2964 0 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%gnorm)
2965 : END IF
2966 : END IF
2967 460 : CALL ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp)
2968 :
2969 : ! The m+1 slot was scratch space for the initial inverse and convergence residual. Preserve
2970 : ! the current accepted point and physical gradient there for the next L-BFGS update.
2971 1154 : DO ispin = 1, nspin
2972 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(m_history + 1)%matrix, &
2973 694 : qs_ot_env(ispin)%matrix_x)
2974 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_e(m_history + 1)%matrix, &
2975 694 : qs_ot_env(ispin)%matrix_gx)
2976 694 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2977 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x_im(m_history + 1)%matrix, &
2978 694 : qs_ot_env(ispin)%matrix_x_im)
2979 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_e_im(m_history + 1)%matrix, &
2980 694 : qs_ot_env(ispin)%matrix_gx_im)
2981 : END IF
2982 694 : IF (qs_ot_env(1)%settings%do_rotation) THEN
2983 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x(m_history + 1)%matrix, &
2984 566 : qs_ot_env(ispin)%rot_mat_x)
2985 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(m_history + 1)%matrix, &
2986 566 : qs_ot_env(ispin)%rot_mat_gx)
2987 566 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
2988 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x_im(m_history + 1)%matrix, &
2989 566 : qs_ot_env(ispin)%rot_mat_x_im)
2990 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e_im(m_history + 1)%matrix, &
2991 566 : qs_ot_env(ispin)%rot_mat_gx_im)
2992 : END IF
2993 : END IF
2994 1154 : IF (do_ener) THEN
2995 2612 : qs_ot_env(ispin)%ener_h_x(m_history + 1, :) = qs_ot_env(ispin)%ener_x
2996 2612 : qs_ot_env(ispin)%ener_h_e(m_history + 1, :) = qs_ot_env(ispin)%ener_gx
2997 : END IF
2998 : END DO
2999 :
3000 : ! Second loop, oldest to newest: z <- V H0 V^T g.
3001 1966 : DO i = nhistory, 1, -1
3002 1506 : beta = 0.0_dp
3003 3754 : DO ispin = 1, nspin
3004 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e(history_index(i))%matrix, &
3005 2248 : qs_ot_env(ispin)%matrix_dx, tmp)
3006 2248 : beta = beta + tmp
3007 2248 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3008 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e_im(history_index(i))%matrix, &
3009 2248 : qs_ot_env(ispin)%matrix_dx_im, tmp)
3010 2248 : beta = beta + tmp
3011 : END IF
3012 2248 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3013 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e(history_index(i))%matrix, &
3014 2152 : qs_ot_env(ispin)%rot_mat_dx, tmp)
3015 2152 : beta = beta + 0.5_dp*tmp
3016 2152 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3017 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e_im(history_index(i))%matrix, &
3018 2152 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
3019 2152 : beta = beta + 0.5_dp*tmp
3020 : END IF
3021 : END IF
3022 6002 : IF (do_ener) THEN
3023 : beta = beta + DOT_PRODUCT(qs_ot_env(ispin)%ener_h_e(history_index(i), :), &
3024 10278 : qs_ot_env(ispin)%ener_dx)
3025 : END IF
3026 : END DO
3027 1506 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, beta)
3028 1506 : beta = rho(i)*beta
3029 4214 : DO ispin = 1, nspin
3030 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, &
3031 : qs_ot_env(ispin)%matrix_h_x(history_index(i))%matrix, &
3032 2248 : alpha_scalar=1.0_dp, beta_scalar=alpha(i) - beta)
3033 2248 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3034 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx_im, &
3035 : qs_ot_env(ispin)%matrix_h_x_im(history_index(i))%matrix, &
3036 2248 : alpha_scalar=1.0_dp, beta_scalar=alpha(i) - beta)
3037 : END IF
3038 2248 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3039 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
3040 : qs_ot_env(ispin)%rot_mat_h_x(history_index(i))%matrix, &
3041 2152 : alpha_scalar=1.0_dp, beta_scalar=alpha(i) - beta)
3042 2152 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3043 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx_im, &
3044 : qs_ot_env(ispin)%rot_mat_h_x_im(history_index(i))%matrix, &
3045 2152 : alpha_scalar=1.0_dp, beta_scalar=alpha(i) - beta)
3046 : END IF
3047 : END IF
3048 3754 : IF (do_ener) THEN
3049 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_dx + &
3050 : (alpha(i) - beta)* &
3051 10278 : qs_ot_env(ispin)%ener_h_x(history_index(i), :)
3052 : END IF
3053 : END DO
3054 : END DO
3055 :
3056 460 : g_dot_z = 0.0_dp
3057 1154 : DO ispin = 1, nspin
3058 694 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
3059 694 : g_dot_z = g_dot_z + tmp
3060 694 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3061 694 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_dx_im, tmp)
3062 694 : g_dot_z = g_dot_z + tmp
3063 : END IF
3064 694 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3065 566 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
3066 566 : g_dot_z = g_dot_z + 0.5_dp*tmp
3067 566 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3068 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
3069 566 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
3070 566 : g_dot_z = g_dot_z + 0.5_dp*tmp
3071 : END IF
3072 : END IF
3073 1848 : IF (do_ener) THEN
3074 : g_dot_z = g_dot_z + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
3075 2612 : qs_ot_env(ispin)%ener_dx)
3076 : END IF
3077 : END DO
3078 460 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_z)
3079 :
3080 : ! Preserve a guaranteed descent direction even if the supplied initial preconditioner is
3081 : ! indefinite or numerical noise corrupts the secant recursion.
3082 460 : IF (.NOT. ieee_is_finite(g_dot_z) .OR. g_dot_z <= 0.0_dp) THEN
3083 0 : qs_ot_env(1)%diis_iter = 1
3084 0 : qs_ot_env(1)%OT_METHOD_FULL = "OT L-SD"
3085 0 : DO ispin = 1, nspin
3086 : CALL ot_apply_channel_preconditioner(qs_ot_env(ispin), &
3087 : qs_ot_env(ispin)%matrix_gx, &
3088 : qs_ot_env(ispin)%matrix_dx, &
3089 : qs_ot_env(ispin)%matrix_gx_im, &
3090 0 : qs_ot_env(ispin)%matrix_dx_im)
3091 0 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3092 0 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
3093 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3094 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, &
3095 0 : qs_ot_env(ispin)%rot_mat_gx_im)
3096 : END IF
3097 : END IF
3098 0 : IF (do_ener) qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_gx
3099 : END DO
3100 0 : g_dot_z = 0.0_dp
3101 0 : DO ispin = 1, nspin
3102 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
3103 0 : g_dot_z = g_dot_z + tmp
3104 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3105 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, qs_ot_env(ispin)%matrix_dx_im, tmp)
3106 0 : g_dot_z = g_dot_z + tmp
3107 : END IF
3108 0 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3109 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
3110 0 : g_dot_z = g_dot_z + 0.5_dp*tmp
3111 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3112 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
3113 0 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
3114 0 : g_dot_z = g_dot_z + 0.5_dp*tmp
3115 : END IF
3116 : END IF
3117 0 : IF (do_ener) THEN
3118 : g_dot_z = g_dot_z + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
3119 0 : qs_ot_env(ispin)%ener_dx)
3120 : END IF
3121 : END DO
3122 0 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_z)
3123 0 : IF (.NOT. ieee_is_finite(g_dot_z) .OR. g_dot_z <= 0.0_dp) THEN
3124 0 : g_dot_z = 0.0_dp
3125 0 : DO ispin = 1, nspin
3126 0 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx)
3127 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3128 0 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_dx_im, qs_ot_env(ispin)%matrix_gx_im)
3129 : END IF
3130 0 : IF (qs_ot_env(ispin)%settings%ot_algorithm == 'TOD') THEN
3131 : ! The P-projected gradient needs an identity-metric projection for this step.
3132 : CALL qs_ot_project_gradient(qs_ot_env(ispin), qs_ot_env(ispin)%matrix_dx, &
3133 0 : qs_ot_env(ispin)%matrix_dx_im, use_preconditioner=.FALSE.)
3134 : END IF
3135 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
3136 0 : g_dot_z = g_dot_z + tmp
3137 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3138 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
3139 0 : qs_ot_env(ispin)%matrix_dx_im, tmp)
3140 0 : g_dot_z = g_dot_z + tmp
3141 : END IF
3142 0 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3143 0 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
3144 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
3145 0 : g_dot_z = g_dot_z + 0.5_dp*tmp
3146 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3147 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx_im, &
3148 0 : qs_ot_env(ispin)%rot_mat_gx_im)
3149 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
3150 0 : qs_ot_env(ispin)%rot_mat_dx_im, tmp)
3151 0 : g_dot_z = g_dot_z + 0.5_dp*tmp
3152 : END IF
3153 : END IF
3154 0 : IF (do_ener) THEN
3155 0 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_gx
3156 : g_dot_z = g_dot_z + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
3157 0 : qs_ot_env(ispin)%ener_dx)
3158 : END IF
3159 : END DO
3160 0 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, g_dot_z)
3161 : END IF
3162 : END IF
3163 :
3164 1154 : DO ispin = 1, nspin
3165 694 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
3166 694 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3167 694 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx_im, -1.0_dp)
3168 : END IF
3169 694 : IF (qs_ot_env(1)%settings%do_rotation) THEN
3170 566 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
3171 566 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3172 566 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx_im, -1.0_dp)
3173 : END IF
3174 : END IF
3175 3498 : IF (do_ener) qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_dx
3176 : END DO
3177 460 : qs_ot_env(1)%gnorm = g_dot_z
3178 460 : qs_ot_env(1)%gradient = -g_dot_z
3179 :
3180 460 : IF (nvariables <= 0.0_dp) qs_ot_env(1)%gradient = 0.0_dp
3181 :
3182 460 : test_down = -g_dot_z
3183 : CALL ot_try_mermin_response_direction( &
3184 460 : qs_ot_env, para_env_inter_kp, qs_ot_env(1)%delta, test_down, use_response_candidate)
3185 460 : IF (use_response_candidate) THEN
3186 : ! A sparse finite-response probe changes the inverse model outside the accumulated
3187 : ! L-BFGS secant space. Restart that space while retaining the accepted point used by the
3188 : ! next update; this avoids mixing a calibrated physical candidate with stale curvature.
3189 8 : qs_ot_env(1)%diis_iter = 1
3190 8 : qs_ot_env(1)%OT_METHOD_FULL = "OT L-R"
3191 8 : qs_ot_env(1)%gnorm = -test_down
3192 8 : qs_ot_env(1)%gradient = test_down
3193 : END IF
3194 :
3195 460 : DEALLOCATE (alpha, rho, history_index)
3196 460 : CALL timestop(handle)
3197 :
3198 1380 : END SUBROUTINE ot_new_lbfgs_direction
3199 :
3200 : ! **************************************************************************************************
3201 : !> \brief Build a preconditioned product residual for history-based OT minimizers.
3202 : !> \param qs_ot_env OT environments for all local spin/k-point channels
3203 : !> \param history_index destination slot in the residual history
3204 : !> \param step_scale scale applied after evaluating the descent product
3205 : !> \param para_env_inter_kp communicator between distributed k-point groups
3206 : !> \param gnorm product of the physical Mermin gradient and preconditioned residual
3207 : ! **************************************************************************************************
3208 51328 : SUBROUTINE ot_build_history_residual(qs_ot_env, history_index, step_scale, &
3209 : para_env_inter_kp, gnorm)
3210 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
3211 : INTEGER, INTENT(IN) :: history_index
3212 : REAL(KIND=dp), INTENT(IN) :: step_scale
3213 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
3214 : REAL(KIND=dp), INTENT(OUT) :: gnorm
3215 :
3216 : LOGICAL :: use_occupation_response
3217 : TYPE(cp_logger_type), POINTER :: logger
3218 :
3219 51328 : use_occupation_response = qs_ot_env(1)%settings%occupation_preconditioner
3220 : CALL ot_build_history_residual_once(qs_ot_env, history_index, step_scale, &
3221 51328 : use_occupation_response, gnorm)
3222 51328 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, gnorm)
3223 :
3224 : ! The occupation response is a coupled orbital/energy approximation. Falling back one block
3225 : ! at a time would change its Schur balance, so replace the complete product residual instead.
3226 51328 : IF (use_occupation_response .AND. gnorm <= 0.0_dp) THEN
3227 2 : CALL ot_build_history_residual_once(qs_ot_env, history_index, step_scale, .FALSE., gnorm)
3228 2 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, gnorm)
3229 : END IF
3230 51328 : IF (gnorm < 0.0_dp) THEN
3231 0 : logger => cp_get_default_logger()
3232 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
3233 : END IF
3234 51328 : END SUBROUTINE ot_build_history_residual
3235 :
3236 : ! **************************************************************************************************
3237 : !> \brief Build one physical or occupation-response product residual.
3238 : !> \param qs_ot_env OT environments for all local spin/k-point channels
3239 : !> \param history_index destination slot in the residual history
3240 : !> \param step_scale scale applied to the stored residual
3241 : !> \param use_occupation_response use the fixed-N occupation-response approximation
3242 : !> \param gnorm product of the physical Mermin gradient and unscaled residual
3243 : !> \param use_preconditioner apply the orbital preconditioner; defaults to true
3244 : ! **************************************************************************************************
3245 51330 : SUBROUTINE ot_build_history_residual_once(qs_ot_env, history_index, step_scale, &
3246 : use_occupation_response, gnorm, use_preconditioner)
3247 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
3248 : INTEGER, INTENT(IN) :: history_index
3249 : REAL(KIND=dp), INTENT(IN) :: step_scale
3250 : LOGICAL, INTENT(IN) :: use_occupation_response
3251 : REAL(KIND=dp), INTENT(OUT) :: gnorm
3252 : LOGICAL, INTENT(IN), OPTIONAL :: use_preconditioner
3253 :
3254 : INTEGER :: ispin
3255 : LOGICAL :: do_ener, do_ks, precondition
3256 : REAL(KIND=dp) :: tmp
3257 : TYPE(dbcsr_type), POINTER :: history_residual_im, orbital_residual, &
3258 : orbital_residual_im
3259 226614 : TYPE(dbcsr_p_type) :: matrix_out(2, SIZE(qs_ot_env))
3260 :
3261 51330 : do_ks = qs_ot_env(1)%settings%ks
3262 51330 : do_ener = qs_ot_env(1)%settings%do_ener
3263 51330 : precondition = .TRUE.
3264 51330 : IF (PRESENT(use_preconditioner)) precondition = use_preconditioner
3265 51330 : gnorm = 0.0_dp
3266 :
3267 51330 : IF (do_ks) THEN
3268 51330 : IF (precondition) THEN
3269 109758 : DO ispin = 1, SIZE(qs_ot_env)
3270 58428 : matrix_out(1, ispin)%matrix => qs_ot_env(ispin)%matrix_h_e(history_index)%matrix
3271 109758 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3272 2752 : matrix_out(2, ispin)%matrix => qs_ot_env(ispin)%matrix_h_e_im(history_index)%matrix
3273 : END IF
3274 : END DO
3275 51330 : CALL ot_precondition_gradients(qs_ot_env, matrix_out, use_occupation_response)
3276 : END IF
3277 109758 : DO ispin = 1, SIZE(qs_ot_env)
3278 58428 : NULLIFY (history_residual_im, orbital_residual_im)
3279 58428 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3280 2752 : history_residual_im => qs_ot_env(ispin)%matrix_h_e_im(history_index)%matrix
3281 : END IF
3282 58428 : IF (use_occupation_response) THEN
3283 148 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%matrix_preconditioned_gx))
3284 148 : orbital_residual => qs_ot_env(ispin)%matrix_preconditioned_gx
3285 148 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3286 148 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%matrix_preconditioned_gx_im))
3287 148 : orbital_residual_im => qs_ot_env(ispin)%matrix_preconditioned_gx_im
3288 : END IF
3289 : ELSE
3290 58280 : orbital_residual => qs_ot_env(ispin)%matrix_gx
3291 58280 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3292 2604 : orbital_residual_im => qs_ot_env(ispin)%matrix_gx_im
3293 : END IF
3294 : END IF
3295 :
3296 58428 : IF (.NOT. precondition) THEN
3297 0 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_e(history_index)%matrix, orbital_residual)
3298 0 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3299 0 : CALL dbcsr_copy(history_residual_im, orbital_residual_im)
3300 : END IF
3301 0 : IF (qs_ot_env(ispin)%settings%ot_algorithm == 'TOD') THEN
3302 : CALL qs_ot_project_gradient(qs_ot_env(ispin), &
3303 : qs_ot_env(ispin)%matrix_h_e(history_index)%matrix, &
3304 0 : history_residual_im, use_preconditioner=.FALSE.)
3305 : END IF
3306 : END IF
3307 :
3308 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
3309 58428 : qs_ot_env(ispin)%matrix_h_e(history_index)%matrix, tmp)
3310 58428 : gnorm = gnorm + tmp
3311 58428 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(history_index)%matrix, step_scale)
3312 58428 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3313 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx_im, &
3314 2752 : qs_ot_env(ispin)%matrix_h_e_im(history_index)%matrix, tmp)
3315 2752 : gnorm = gnorm + tmp
3316 2752 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e_im(history_index)%matrix, step_scale)
3317 : END IF
3318 :
3319 168186 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
3320 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(history_index)%matrix, &
3321 750 : qs_ot_env(ispin)%rot_mat_gx)
3322 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
3323 750 : qs_ot_env(ispin)%rot_mat_h_e(history_index)%matrix, tmp)
3324 750 : gnorm = gnorm + 0.5_dp*tmp
3325 750 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e(history_index)%matrix, step_scale)
3326 750 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3327 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e_im(history_index)%matrix, &
3328 392 : qs_ot_env(ispin)%rot_mat_gx_im)
3329 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx_im, &
3330 392 : qs_ot_env(ispin)%rot_mat_h_e_im(history_index)%matrix, tmp)
3331 392 : gnorm = gnorm + 0.5_dp*tmp
3332 392 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e_im(history_index)%matrix, step_scale)
3333 : END IF
3334 : END IF
3335 : END DO
3336 : END IF
3337 :
3338 51330 : IF (do_ener) THEN
3339 506 : DO ispin = 1, SIZE(qs_ot_env)
3340 328 : IF (use_occupation_response) THEN
3341 148 : CPASSERT(ASSOCIATED(qs_ot_env(ispin)%ener_preconditioned_gx))
3342 : qs_ot_env(ispin)%ener_h_e(history_index, :) = &
3343 1652 : qs_ot_env(ispin)%ener_preconditioned_gx
3344 : ELSE
3345 1424 : qs_ot_env(ispin)%ener_h_e(history_index, :) = qs_ot_env(ispin)%ener_gx
3346 : END IF
3347 : gnorm = gnorm + DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, &
3348 3076 : qs_ot_env(ispin)%ener_h_e(history_index, :))
3349 : qs_ot_env(ispin)%ener_h_e(history_index, :) = &
3350 3254 : step_scale*qs_ot_env(ispin)%ener_h_e(history_index, :)
3351 : END DO
3352 : END IF
3353 51330 : END SUBROUTINE ot_build_history_residual_once
3354 :
3355 : ! **************************************************************************************************
3356 : !> \brief ...
3357 : !> \param qs_ot_env ...
3358 : !> \param para_env_inter_kp communicator between distributed k-point groups
3359 : ! **************************************************************************************************
3360 51136 : SUBROUTINE ot_diis_step(qs_ot_env, para_env_inter_kp)
3361 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
3362 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
3363 :
3364 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_diis_step'
3365 :
3366 : INTEGER :: component, diis_bound, diis_m, handle, i, info, &
3367 : ispin, j, m, nmat, nspin
3368 : LOGICAL :: do_ener, do_ks, do_ot_sd, local_history
3369 : REAL(KIND=dp) :: nvariables, overlap, tmp, tr_xnew_gx, &
3370 : tr_xold_gx
3371 51136 : REAL(dp), ALLOCATABLE :: coefficients(:), dots(:), dots_new(:), weights(:)
3372 51136 : TYPE(dbcsr_p_type), ALLOCATABLE :: positions(:), gradients(:), history(:, :)
3373 : TYPE(cp_logger_type), POINTER :: logger
3374 : TYPE(mp_comm_type) :: group
3375 :
3376 51136 : CALL timeset(routineN, handle)
3377 :
3378 51136 : logger => cp_get_default_logger()
3379 :
3380 51136 : do_ks = qs_ot_env(1)%settings%ks
3381 51136 : do_ener = qs_ot_env(1)%settings%do_ener
3382 51136 : nspin = SIZE(qs_ot_env)
3383 :
3384 51136 : diis_m = qs_ot_env(1)%settings%diis_m
3385 :
3386 51136 : IF (qs_ot_env(1)%diis_iter < diis_m) THEN
3387 30882 : diis_bound = qs_ot_env(1)%diis_iter + 1
3388 : ELSE
3389 20254 : diis_bound = diis_m
3390 : END IF
3391 :
3392 51136 : j = MOD(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
3393 :
3394 : ! Keep component order identical to the scalar product: real, imaginary, then rotation.
3395 51136 : nmat = 0
3396 51136 : local_history = do_ks
3397 51136 : IF (do_ks) THEN
3398 109342 : DO ispin = 1, nspin
3399 58206 : m = 1
3400 58206 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) m = 2
3401 58206 : nmat = nmat + m
3402 58206 : IF (qs_ot_env(ispin)%settings%do_rotation) nmat = nmat + m
3403 : local_history = local_history .AND. qs_ot_env(ispin)%has_complex_kpoint_state .AND. &
3404 165650 : .NOT. qs_ot_env(ispin)%settings%do_rotation .AND. qs_ot_env(ispin)%para_env%num_pe == 1
3405 : END DO
3406 : END IF
3407 0 : ALLOCATE (positions(nmat), gradients(nmat), history(nmat, 2*diis_bound), &
3408 1850628 : weights(nmat), dots(nmat), dots_new(nmat), coefficients(2*diis_bound))
3409 51136 : component = 0
3410 51136 : IF (do_ks) THEN
3411 109342 : DO ispin = 1, nspin
3412 : CALL add_component(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_gx, &
3413 58206 : qs_ot_env(ispin)%matrix_h_e, qs_ot_env(ispin)%matrix_h_x, 1.0_dp)
3414 58206 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3415 : CALL add_component(qs_ot_env(ispin)%matrix_x_im, qs_ot_env(ispin)%matrix_gx_im, &
3416 2696 : qs_ot_env(ispin)%matrix_h_e_im, qs_ot_env(ispin)%matrix_h_x_im, 1.0_dp)
3417 : END IF
3418 109342 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
3419 : CALL add_component(qs_ot_env(ispin)%rot_mat_x, qs_ot_env(ispin)%rot_mat_gx, &
3420 694 : qs_ot_env(ispin)%rot_mat_h_e, qs_ot_env(ispin)%rot_mat_h_x, 0.5_dp)
3421 694 : IF (qs_ot_env(ispin)%has_complex_kpoint_state) THEN
3422 : CALL add_component(qs_ot_env(ispin)%rot_mat_x_im, qs_ot_env(ispin)%rot_mat_gx_im, &
3423 336 : qs_ot_env(ispin)%rot_mat_h_e_im, qs_ot_env(ispin)%rot_mat_h_x_im, 0.5_dp)
3424 : END IF
3425 : END IF
3426 : END DO
3427 : END IF
3428 :
3429 : ! copy the position and the error vector in the diis buffers
3430 :
3431 51136 : IF (nmat > 0) CALL dbcsr_get_info(gradients(1)%matrix, group=group)
3432 113068 : DO m = 1, nmat
3433 113068 : CALL dbcsr_copy(history(m, diis_bound + j)%matrix, positions(m)%matrix)
3434 : END DO
3435 51136 : IF (do_ener) THEN
3436 422 : DO ispin = 1, nspin
3437 3002 : qs_ot_env(ispin)%ener_h_x(j, :) = qs_ot_env(ispin)%ener_x(:)
3438 : END DO
3439 : END IF
3440 : CALL ot_build_history_residual(qs_ot_env, j, -qs_ot_env(1)%ds_min, &
3441 51136 : para_env_inter_kp, qs_ot_env(1)%gnorm)
3442 : CALL ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp, &
3443 51136 : residual_product=qs_ot_env(1)%gnorm)
3444 : ! Handling the case of no free variables to optimize
3445 51136 : IF (nvariables > 0.0_dp) THEN
3446 51110 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
3447 : ELSE
3448 26 : qs_ot_env(1)%gradient = 0.0_dp
3449 : END IF
3450 :
3451 : ! make the diis matrix and solve it
3452 304641 : DO i = 1, diis_bound
3453 : ! Form the DIIS metric from the gradient and stored steps.
3454 253505 : qs_ot_env(1)%ls_diis(i, j) = 0.0_dp
3455 253505 : CALL dbcsr_dot_batch(gradients, history(:, i), dots, comm=group)
3456 556566 : DO m = 1, nmat
3457 556566 : qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*weights(m)*dots(m)
3458 : END DO
3459 253505 : IF (do_ener) THEN
3460 2046 : DO ispin = 1, nspin
3461 14166 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_h_e(i, :))
3462 2046 : qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*tmp
3463 : END DO
3464 : END IF
3465 253505 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, qs_ot_env(1)%ls_diis(i, j))
3466 253505 : qs_ot_env(1)%ls_diis(j, i) = qs_ot_env(1)%ls_diis(i, j)
3467 253505 : qs_ot_env(1)%ls_diis(i, diis_bound + 1) = 1.0_dp
3468 253505 : qs_ot_env(1)%ls_diis(diis_bound + 1, i) = 1.0_dp
3469 304641 : qs_ot_env(1)%c_diis(i) = 0.0_dp
3470 : END DO
3471 51136 : qs_ot_env(1)%ls_diis(diis_bound + 1, diis_bound + 1) = 0.0_dp
3472 51136 : qs_ot_env(1)%c_diis(diis_bound + 1) = 1.0_dp
3473 : ! put in buffer, dgesv destroys
3474 3769172 : qs_ot_env(1)%lss_diis = qs_ot_env(1)%ls_diis
3475 :
3476 : CALL DGESV(diis_bound + 1, 1, qs_ot_env(1)%lss_diis, diis_m + 1, qs_ot_env(1)%ipivot, &
3477 51136 : qs_ot_env(1)%c_diis, diis_m + 1, info)
3478 :
3479 51136 : IF (info /= 0) THEN
3480 8 : do_ot_sd = .TRUE.
3481 8 : WRITE (cp_logger_get_default_unit_nr(logger), *) "Singular DIIS matrix"
3482 : ELSE
3483 51128 : do_ot_sd = .FALSE.
3484 304593 : coefficients(1:diis_bound) = qs_ot_env(1)%c_diis(1:diis_bound)
3485 304593 : coefficients(diis_bound + 1:2*diis_bound) = qs_ot_env(1)%c_diis(1:diis_bound)
3486 51128 : IF (local_history) THEN
3487 772 : CALL dbcsr_sum_many(positions, history, coefficients)
3488 : ELSE
3489 108468 : DO m = 1, nmat
3490 58112 : CALL dbcsr_set(positions(m)%matrix, 0.0_dp)
3491 681054 : DO i = 1, 2*diis_bound
3492 630698 : CALL dbcsr_add(positions(m)%matrix, history(m, i)%matrix, 1.0_dp, coefficients(i))
3493 : END DO
3494 : END DO
3495 : END IF
3496 51128 : IF (do_ener) THEN
3497 422 : DO ispin = 1, nspin
3498 2852 : qs_ot_env(ispin)%ener_x(:) = 0.0_dp
3499 1598 : DO i = 1, diis_bound
3500 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
3501 14438 : + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_e(i, :)
3502 : END DO
3503 1748 : DO i = 1, diis_bound
3504 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
3505 14438 : + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_x(i, :)
3506 : END DO
3507 : END DO
3508 : END IF
3509 51128 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
3510 51128 : IF (qs_ot_env(1)%settings%safer_diis) THEN
3511 : ! now, final check, is the step in fact in the direction of the -gradient ?
3512 : ! if not we're walking towards a sadle point, and should avoid that
3513 : ! the direction of the step is x_new-x_old
3514 51128 : tr_xold_gx = 0.0_dp
3515 51128 : tr_xnew_gx = 0.0_dp
3516 51128 : CALL dbcsr_dot_batch(history(:, diis_bound + j), gradients, dots, comm=group)
3517 51128 : CALL dbcsr_dot_batch(positions, gradients, dots_new, comm=group)
3518 113036 : DO m = 1, nmat
3519 61908 : tr_xold_gx = tr_xold_gx + weights(m)*dots(m)
3520 113036 : tr_xnew_gx = tr_xnew_gx + weights(m)*dots_new(m)
3521 : END DO
3522 51128 : IF (do_ener) THEN
3523 422 : DO ispin = 1, nspin
3524 2852 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_x(j, :), qs_ot_env(ispin)%ener_gx(:))
3525 272 : tr_xold_gx = tr_xold_gx + tmp
3526 2852 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_x(:), qs_ot_env(ispin)%ener_gx(:))
3527 422 : tr_xnew_gx = tr_xnew_gx + tmp
3528 : END DO
3529 : END IF
3530 51128 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, tr_xold_gx)
3531 51128 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, tr_xnew_gx)
3532 51128 : overlap = (tr_xnew_gx - tr_xold_gx)
3533 : ! OK, bad luck, take a SD step along the preconditioned gradient
3534 51128 : IF (overlap > 0.0_dp) THEN
3535 : do_ot_sd = .TRUE.
3536 : END IF
3537 : END IF
3538 : END IF
3539 :
3540 : IF (do_ot_sd) THEN
3541 1118 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
3542 2430 : DO m = 1, nmat
3543 1312 : CALL dbcsr_set(positions(m)%matrix, 0.0_dp)
3544 1312 : CALL dbcsr_add(positions(m)%matrix, history(m, j)%matrix, 1.0_dp, 1.0_dp)
3545 2430 : CALL dbcsr_add(positions(m)%matrix, history(m, diis_bound + j)%matrix, 1.0_dp, 1.0_dp)
3546 : END DO
3547 1118 : IF (do_ener) THEN
3548 12 : DO ispin = 1, nspin
3549 184 : qs_ot_env(ispin)%ener_x(:) = 0._dp
3550 184 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_e(j, :)
3551 188 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_x(j, :)
3552 : END DO
3553 : END IF
3554 : END IF
3555 :
3556 153408 : CALL timestop(handle)
3557 :
3558 : CONTAINS
3559 :
3560 : ! **************************************************************************************************
3561 : !> \brief Collect one DIIS component without copying matrix data.
3562 : !> \param x current position
3563 : !> \param gx current gradient
3564 : !> \param he error history
3565 : !> \param hx position history
3566 : !> \param weight metric weight of this component
3567 : ! **************************************************************************************************
3568 61932 : SUBROUTINE add_component(x, gx, he, hx, weight)
3569 : TYPE(dbcsr_type), POINTER :: x, gx
3570 : TYPE(dbcsr_p_type), INTENT(IN) :: he(:), hx(:)
3571 : REAL(dp), INTENT(IN) :: weight
3572 :
3573 61932 : component = component + 1
3574 61932 : positions(component)%matrix => x
3575 61932 : gradients(component)%matrix => gx
3576 364993 : history(component, 1:diis_bound) = he(1:diis_bound)
3577 364993 : history(component, diis_bound + 1:2*diis_bound) = hx(1:diis_bound)
3578 61932 : weights(component) = weight
3579 61932 : END SUBROUTINE add_component
3580 :
3581 : END SUBROUTINE ot_diis_step
3582 :
3583 : ! **************************************************************************************************
3584 : !> \brief Energy minimizer by Broyden's method
3585 : !> \param qs_ot_env variable to control minimizer behaviour
3586 : !> \param para_env_inter_kp communicator between distributed k-point groups
3587 : !> \author Kurt Baarman (09.2010)
3588 : ! **************************************************************************************************
3589 192 : SUBROUTINE ot_broyden_step(qs_ot_env, para_env_inter_kp)
3590 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
3591 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
3592 :
3593 : INTEGER, PARAMETER :: broyden_gradient = 4, &
3594 : broyden_position = 1, &
3595 : broyden_random = 3, &
3596 : broyden_residual = 2
3597 : INTEGER :: diis_bound, diis_m, i, j, k, nspin
3598 192 : INTEGER, ALLOCATABLE, DIMENSION(:) :: circ_index
3599 : LOGICAL :: adaptive_sigma, do_ener, do_ks, &
3600 : do_rotation, enable_flip, forget_history
3601 : REAL(KIND=dp) :: beta, eta, gamma, nvariables, omega, &
3602 : sigma, sigma_dec, sigma_min, tmp, tmp2
3603 192 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f, x
3604 192 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: G, S
3605 192 : eta = qs_ot_env(1)%settings%broyden_eta
3606 192 : omega = qs_ot_env(1)%settings%broyden_omega
3607 192 : sigma_dec = qs_ot_env(1)%settings%broyden_sigma_decrease
3608 192 : sigma_min = qs_ot_env(1)%settings%broyden_sigma_min
3609 192 : forget_history = qs_ot_env(1)%settings%broyden_forget_history
3610 192 : adaptive_sigma = qs_ot_env(1)%settings%broyden_adaptive_sigma
3611 192 : enable_flip = qs_ot_env(1)%settings%broyden_enable_flip
3612 192 : do_ks = qs_ot_env(1)%settings%ks
3613 192 : do_ener = qs_ot_env(1)%settings%do_ener
3614 192 : do_rotation = qs_ot_env(1)%settings%do_rotation
3615 :
3616 192 : beta = qs_ot_env(1)%settings%broyden_beta
3617 192 : gamma = qs_ot_env(1)%settings%broyden_gamma
3618 192 : IF (adaptive_sigma) THEN
3619 192 : IF (qs_ot_env(1)%broyden_adaptive_sigma < 0.0_dp) THEN
3620 6 : sigma = qs_ot_env(1)%settings%broyden_sigma
3621 : ELSE
3622 : sigma = qs_ot_env(1)%broyden_adaptive_sigma
3623 : END IF
3624 : ELSE
3625 0 : sigma = qs_ot_env(1)%settings%broyden_sigma
3626 : END IF
3627 :
3628 192 : IF (.NOT. do_ks) CPABORT("BROYDEN currently requires OT orbital variables")
3629 :
3630 192 : nspin = SIZE(qs_ot_env)
3631 :
3632 192 : diis_m = qs_ot_env(1)%settings%diis_m
3633 :
3634 192 : IF (qs_ot_env(1)%diis_iter < diis_m) THEN
3635 90 : diis_bound = qs_ot_env(1)%diis_iter + 1
3636 : ELSE
3637 102 : diis_bound = diis_m
3638 : END IF
3639 :
3640 : ! We want x:s, f:s and one random vector
3641 192 : k = 2*diis_bound + 1
3642 768 : ALLOCATE (S(k, k))
3643 576 : ALLOCATE (G(k, k))
3644 576 : ALLOCATE (f(k))
3645 384 : ALLOCATE (x(k))
3646 576 : ALLOCATE (circ_index(diis_bound))
3647 192 : G = 0.0_dp
3648 2500 : DO i = 1, k
3649 2500 : G(i, i) = sigma
3650 : END DO
3651 192 : S = 0.0_dp
3652 :
3653 192 : j = MOD(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
3654 :
3655 192 : CALL broyden_copy_current_to_history(j)
3656 : CALL ot_build_history_residual(qs_ot_env, j, -1.0_dp, para_env_inter_kp, &
3657 192 : qs_ot_env(1)%gnorm)
3658 192 : IF (qs_ot_env(1)%settings%occupation_preconditioner) THEN
3659 26 : CALL broyden_bound_occupation_response(j, qs_ot_env(1)%gnorm)
3660 : END IF
3661 :
3662 : CALL ot_update_delta(qs_ot_env, nvariables, para_env_inter_kp, &
3663 192 : residual_product=qs_ot_env(1)%gnorm)
3664 :
3665 : ! Handling the case of no free variables to optimize
3666 192 : IF (nvariables > 0.0_dp) THEN
3667 192 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
3668 : ELSE
3669 0 : qs_ot_env(1)%gradient = 0.0_dp
3670 : END IF
3671 :
3672 192 : IF (diis_bound == diis_m) THEN
3673 912 : DO i = 1, diis_bound
3674 912 : circ_index(i) = MOD(j + i - 1, diis_m) + 1
3675 : END DO
3676 : ELSE
3677 338 : DO i = 1, diis_bound
3678 338 : circ_index(i) = i
3679 : END DO
3680 : END IF
3681 :
3682 192 : S = 0.0_dp
3683 192 : CALL broyden_randomize_current()
3684 1250 : DO i = 1, diis_bound
3685 : CALL broyden_product_dot(broyden_position, circ_index(i), &
3686 1058 : broyden_position, circ_index(i), S(i, i))
3687 : CALL broyden_product_dot(broyden_residual, circ_index(i), &
3688 : broyden_residual, circ_index(i), &
3689 1058 : S(i + diis_bound, i + diis_bound))
3690 : ! Preserve the established Broyden model, whose auxiliary vector is orthogonal to history.
3691 1058 : S(i, 2*diis_bound + 1) = 0.0_dp
3692 1058 : S(2*diis_bound + 1, i) = 0.0_dp
3693 1058 : S(i + diis_bound, 2*diis_bound + 1) = 0.0_dp
3694 1058 : S(2*diis_bound + 1, i + diis_bound) = 0.0_dp
3695 3874 : DO k = i + 1, diis_bound
3696 : CALL broyden_product_dot(broyden_position, circ_index(i), &
3697 2816 : broyden_position, circ_index(k), S(i, k))
3698 2816 : S(k, i) = S(i, k)
3699 : CALL broyden_product_dot(broyden_residual, circ_index(i), &
3700 : broyden_residual, circ_index(k), &
3701 2816 : S(diis_bound + i, diis_bound + k))
3702 3874 : S(diis_bound + k, diis_bound + i) = S(diis_bound + i, diis_bound + k)
3703 : END DO
3704 7940 : DO k = 1, diis_bound
3705 : CALL broyden_product_dot(broyden_position, circ_index(i), &
3706 : broyden_residual, circ_index(k), &
3707 6690 : S(i, k + diis_bound))
3708 7748 : S(k + diis_bound, i) = S(i, k + diis_bound)
3709 : END DO
3710 : END DO
3711 : CALL broyden_product_dot(broyden_random, 0, broyden_random, 0, &
3712 192 : S(2*diis_bound + 1, 2*diis_bound + 1))
3713 :
3714 : ! normalize
3715 192 : k = 2*diis_bound + 1
3716 192 : tmp = SQRT(S(k, k))
3717 2500 : S(k, :) = S(k, :)/tmp
3718 2500 : S(:, k) = S(:, k)/tmp
3719 :
3720 192 : IF (diis_bound > 1) THEN
3721 176 : tmp2 = 0.0_dp
3722 176 : i = diis_bound
3723 : CALL broyden_product_dot(broyden_position, circ_index(i), &
3724 176 : broyden_residual, circ_index(i), tmp)
3725 176 : tmp2 = tmp2 + tmp
3726 : CALL broyden_product_dot(broyden_position, circ_index(i - 1), &
3727 176 : broyden_residual, circ_index(i), tmp)
3728 176 : tmp2 = tmp2 - tmp
3729 : CALL broyden_product_dot(broyden_position, circ_index(i), &
3730 176 : broyden_residual, circ_index(i - 1), tmp)
3731 176 : tmp2 = tmp2 - tmp
3732 : CALL broyden_product_dot(broyden_position, circ_index(i - 1), &
3733 176 : broyden_residual, circ_index(i - 1), tmp)
3734 176 : tmp2 = tmp2 + tmp
3735 176 : qs_ot_env(1)%c_broy(i - 1) = tmp2
3736 : END IF
3737 :
3738 192 : qs_ot_env(1)%energy_h(j) = qs_ot_env(1)%etotal
3739 :
3740 : ! If we went uphill, do backtracking line search
3741 1250 : i = MINLOC(qs_ot_env(1)%energy_h(1:diis_bound), dim=1)
3742 192 : IF (i /= j) THEN
3743 26 : sigma = sigma_dec*sigma
3744 26 : qs_ot_env(1)%OT_METHOD_FULL = "OT BTRK"
3745 26 : CALL broyden_set_current_zero()
3746 26 : CALL broyden_add_history(broyden_position, i, 1.0_dp, 1.0_dp - gamma)
3747 26 : CALL broyden_add_history(broyden_position, circ_index(diis_bound), 1.0_dp, gamma)
3748 : ELSE
3749 : ! Construct G
3750 876 : DO i = 2, diis_bound
3751 710 : f = 0.0_dp
3752 710 : x = 0.0_dp
3753 : ! f is df_i
3754 710 : x(i) = 1.0_dp
3755 710 : x(i - 1) = -1.0_dp
3756 : ! x is dx_i
3757 710 : f(diis_bound + i) = 1.0_dp
3758 710 : f(diis_bound + i - 1) = -1.0_dp
3759 710 : tmp = 1.0_dp
3760 : ! We want a pos def Hessian
3761 710 : IF (enable_flip) THEN
3762 710 : IF (qs_ot_env(1)%c_broy(i - 1) > 0) THEN
3763 : !qs_ot_env(1)%OT_METHOD_FULL="OT FLIP"
3764 2 : tmp = -1.0_dp
3765 : END IF
3766 : END IF
3767 :
3768 : ! get dx-Gdf
3769 439280 : x(:) = tmp*x - MATMUL(G, f)
3770 : ! dfSdf
3771 : ! we calculate matmul(S, f) twice. They're small...
3772 439280 : tmp = DOT_PRODUCT(f, MATMUL(S, f))
3773 : ! NOTE THAT S IS SYMMETRIC !!!
3774 439280 : f(:) = MATMUL(S, f)/tmp
3775 : ! the spread is an outer vector product
3776 149856 : G(:, :) = G + SPREAD(x, dim=2, ncopies=SIZE(f))*SPREAD(f, dim=1, ncopies=SIZE(x))
3777 : END DO
3778 166 : f = 0.0_dp
3779 166 : f(2*diis_bound) = 1.0_dp
3780 105338 : x(:) = -beta*MATMUL(G, f)
3781 :
3782 : ! OK, add the vectors now, this sums up to the proposed step
3783 166 : CALL broyden_set_current_zero()
3784 1042 : DO i = 1, diis_bound
3785 : CALL broyden_add_history(broyden_residual, circ_index(i), 1.0_dp, &
3786 1042 : -x(i + diis_bound))
3787 : END DO
3788 1042 : DO i = 1, diis_bound
3789 1042 : CALL broyden_add_history(broyden_position, circ_index(i), 1.0_dp, x(i))
3790 : END DO
3791 :
3792 166 : IF (adaptive_sigma) THEN
3793 166 : tmp = new_sigma(G, S, diis_bound)
3794 : !tmp = tmp * qs_ot_env ( 1 ) % settings % broyden_sigma
3795 166 : tmp = tmp*eta
3796 166 : sigma = MIN(omega*sigma, tmp)
3797 : END IF
3798 :
3799 : ! compute the inner product of direction of the step and gradient
3800 166 : CALL broyden_product_dot(broyden_gradient, 0, broyden_random, 0, tmp)
3801 :
3802 : ! If the direction is not a descent direction, reverse it before adding the current point.
3803 166 : IF (tmp >= 0.0_dp) THEN
3804 8 : qs_ot_env(1)%OT_METHOD_FULL = "OT TURN"
3805 : ! A non-descent product step invalidates the secants that couple orbitals,
3806 : ! rotations, and auxiliary energies. Retain the legacy opt-in behavior
3807 : ! for orbital-only OT, but always restart an inconsistent Mermin history.
3808 8 : IF (broyden_history_restart_required(.TRUE., forget_history, do_ener)) THEN
3809 0 : qs_ot_env(1)%diis_iter = 0
3810 : END IF
3811 8 : sigma = sigma*sigma_dec
3812 8 : CALL broyden_add_history(broyden_position, circ_index(diis_bound), -1.0_dp, 1.0_dp)
3813 : ELSE
3814 158 : CALL broyden_add_history(broyden_position, circ_index(diis_bound), 1.0_dp, 1.0_dp)
3815 : END IF
3816 : END IF
3817 :
3818 : ! get rid of S, G, f, x, circ_index for next round
3819 192 : DEALLOCATE (S, G, f, x, circ_index)
3820 :
3821 : ! update for next round
3822 192 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
3823 384 : qs_ot_env(1)%broyden_adaptive_sigma = MAX(sigma, sigma_min)
3824 :
3825 : CONTAINS
3826 :
3827 : ! **************************************************************************************************
3828 : !> \brief Bound the occupation-response direction against the conventional OT H0 direction.
3829 : !> \param history_index response history slot
3830 : !> \param gnorm physical gradient product with the resulting positive direction
3831 : ! **************************************************************************************************
3832 26 : SUBROUTINE broyden_bound_occupation_response(history_index, gnorm)
3833 : INTEGER, INTENT(IN) :: history_index
3834 : REAL(KIND=dp), INTENT(INOUT) :: gnorm
3835 :
3836 : INTEGER :: ichannel
3837 : LOGICAL :: valid_response
3838 : REAL(KIND=dp) :: g_dot_h0_g, h0_g_norm_sq, p_dot_h0, response_delta_norm_sq, &
3839 : response_norm_sq, response_scale, response_weight
3840 :
3841 26 : CALL broyden_build_h0_current()
3842 26 : CALL broyden_product_dot(broyden_gradient, 0, broyden_random, 0, g_dot_h0_g)
3843 26 : CALL broyden_product_dot(broyden_random, 0, broyden_random, 0, h0_g_norm_sq)
3844 : CALL broyden_product_dot(broyden_residual, history_index, &
3845 26 : broyden_residual, history_index, response_norm_sq)
3846 : CALL broyden_product_dot(broyden_residual, history_index, &
3847 26 : broyden_random, 0, p_dot_h0)
3848 26 : p_dot_h0 = -p_dot_h0
3849 26 : response_delta_norm_sq = response_norm_sq + h0_g_norm_sq - 2.0_dp*p_dot_h0
3850 : valid_response = ieee_is_finite(gnorm) .AND. ieee_is_finite(response_norm_sq) .AND. &
3851 : ieee_is_finite(g_dot_h0_g) .AND. ieee_is_finite(h0_g_norm_sq) .AND. &
3852 : ieee_is_finite(response_delta_norm_sq) .AND. gnorm > 0.0_dp .AND. &
3853 : g_dot_h0_g > 0.0_dp .AND. h0_g_norm_sq > TINY(1.0_dp) .AND. &
3854 26 : response_delta_norm_sq >= 0.0_dp
3855 52 : response_scale = 1.0_dp
3856 26 : response_weight = 0.0_dp
3857 : IF (valid_response) THEN
3858 : ! Keep the response correction inside a product-norm trust ball around H0*g.
3859 : response_weight = MIN(1.0_dp, SQRT(h0_g_norm_sq/ &
3860 26 : MAX(response_delta_norm_sq, TINY(1.0_dp))))
3861 : END IF
3862 : ! Broyden secants require a consistent residual map. Use the state-dependent occupation
3863 : ! response only to seed the first direction, then let the history learn from frozen H0*g.
3864 26 : IF (qs_ot_env(1)%diis_iter > 0) response_weight = 0.0_dp
3865 26 : IF (.NOT. valid_response) THEN
3866 0 : response_scale = 0.0_dp
3867 0 : response_weight = 0.0_dp
3868 : END IF
3869 :
3870 78 : DO ichannel = 1, nspin
3871 : CALL dbcsr_scale(qs_ot_env(ichannel)%matrix_h_e(history_index)%matrix, &
3872 52 : response_weight*response_scale)
3873 : CALL dbcsr_add(qs_ot_env(ichannel)%matrix_h_e(history_index)%matrix, &
3874 52 : qs_ot_env(ichannel)%matrix_x, 1.0_dp, -(1.0_dp - response_weight))
3875 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3876 : CALL dbcsr_scale(qs_ot_env(ichannel)%matrix_h_e_im(history_index)%matrix, &
3877 52 : response_weight*response_scale)
3878 : CALL dbcsr_add(qs_ot_env(ichannel)%matrix_h_e_im(history_index)%matrix, &
3879 : qs_ot_env(ichannel)%matrix_x_im, 1.0_dp, &
3880 52 : -(1.0_dp - response_weight))
3881 : END IF
3882 52 : IF (do_rotation) THEN
3883 : CALL dbcsr_scale(qs_ot_env(ichannel)%rot_mat_h_e(history_index)%matrix, &
3884 52 : response_weight*response_scale)
3885 : CALL dbcsr_add(qs_ot_env(ichannel)%rot_mat_h_e(history_index)%matrix, &
3886 : qs_ot_env(ichannel)%rot_mat_x, 1.0_dp, &
3887 52 : -(1.0_dp - response_weight))
3888 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3889 : CALL dbcsr_scale(qs_ot_env(ichannel)%rot_mat_h_e_im(history_index)%matrix, &
3890 52 : response_weight*response_scale)
3891 : CALL dbcsr_add(qs_ot_env(ichannel)%rot_mat_h_e_im(history_index)%matrix, &
3892 : qs_ot_env(ichannel)%rot_mat_x_im, 1.0_dp, &
3893 52 : -(1.0_dp - response_weight))
3894 : END IF
3895 : END IF
3896 78 : IF (do_ener) THEN
3897 : qs_ot_env(ichannel)%ener_h_e(history_index, :) = &
3898 : response_weight*response_scale*qs_ot_env(ichannel)%ener_h_e(history_index, :) - &
3899 208 : (1.0_dp - response_weight)*qs_ot_env(ichannel)%ener_x
3900 : END IF
3901 : END DO
3902 26 : IF (valid_response) THEN
3903 : gnorm = response_weight*response_scale*gnorm + &
3904 26 : (1.0_dp - response_weight)*g_dot_h0_g
3905 : ELSE
3906 0 : gnorm = g_dot_h0_g
3907 : END IF
3908 26 : END SUBROUTINE broyden_bound_occupation_response
3909 :
3910 : ! **************************************************************************************************
3911 : !> \brief Store the conventional positive OT H0 direction in the current product coordinate.
3912 : ! **************************************************************************************************
3913 26 : SUBROUTINE broyden_build_h0_current()
3914 : INTEGER :: ichannel
3915 :
3916 78 : DO ichannel = 1, nspin
3917 : CALL ot_apply_channel_preconditioner(qs_ot_env(ichannel), &
3918 : qs_ot_env(ichannel)%matrix_gx, &
3919 : qs_ot_env(ichannel)%matrix_x, &
3920 : qs_ot_env(ichannel)%matrix_gx_im, &
3921 52 : qs_ot_env(ichannel)%matrix_x_im)
3922 52 : IF (do_rotation) THEN
3923 52 : CALL dbcsr_copy(qs_ot_env(ichannel)%rot_mat_x, qs_ot_env(ichannel)%rot_mat_gx)
3924 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3925 : CALL dbcsr_copy(qs_ot_env(ichannel)%rot_mat_x_im, &
3926 52 : qs_ot_env(ichannel)%rot_mat_gx_im)
3927 : END IF
3928 : END IF
3929 234 : IF (do_ener) qs_ot_env(ichannel)%ener_x = qs_ot_env(ichannel)%ener_gx
3930 : END DO
3931 26 : END SUBROUTINE broyden_build_h0_current
3932 :
3933 : ! **************************************************************************************************
3934 : !> \brief Copy all current product coordinates into a Broyden history slot.
3935 : !> \param history_index destination history slot
3936 : ! **************************************************************************************************
3937 192 : SUBROUTINE broyden_copy_current_to_history(history_index)
3938 : INTEGER, INTENT(IN) :: history_index
3939 :
3940 : INTEGER :: ichannel
3941 :
3942 410 : DO ichannel = 1, nspin
3943 : CALL dbcsr_copy(qs_ot_env(ichannel)%matrix_h_x(history_index)%matrix, &
3944 218 : qs_ot_env(ichannel)%matrix_x)
3945 218 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3946 : CALL dbcsr_copy(qs_ot_env(ichannel)%matrix_h_x_im(history_index)%matrix, &
3947 52 : qs_ot_env(ichannel)%matrix_x_im)
3948 : END IF
3949 218 : IF (do_rotation) THEN
3950 : CALL dbcsr_copy(qs_ot_env(ichannel)%rot_mat_h_x(history_index)%matrix, &
3951 52 : qs_ot_env(ichannel)%rot_mat_x)
3952 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3953 : CALL dbcsr_copy(qs_ot_env(ichannel)%rot_mat_h_x_im(history_index)%matrix, &
3954 52 : qs_ot_env(ichannel)%rot_mat_x_im)
3955 : END IF
3956 : END IF
3957 410 : IF (do_ener) THEN
3958 208 : qs_ot_env(ichannel)%ener_h_x(history_index, :) = qs_ot_env(ichannel)%ener_x
3959 : END IF
3960 : END DO
3961 192 : END SUBROUTINE broyden_copy_current_to_history
3962 :
3963 : ! **************************************************************************************************
3964 : !> \brief Fill the current product coordinate with a deterministic auxiliary vector.
3965 : ! **************************************************************************************************
3966 192 : SUBROUTINE broyden_randomize_current()
3967 : INTEGER :: ichannel, icoef
3968 :
3969 410 : DO ichannel = 1, nspin
3970 218 : CALL dbcsr_init_random(qs_ot_env(ichannel)%matrix_x)
3971 218 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3972 52 : CALL dbcsr_init_random(qs_ot_env(ichannel)%matrix_x_im)
3973 : END IF
3974 218 : IF (do_rotation) THEN
3975 52 : CALL dbcsr_init_random(qs_ot_env(ichannel)%rot_mat_x)
3976 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3977 52 : CALL dbcsr_init_random(qs_ot_env(ichannel)%rot_mat_x_im)
3978 : END IF
3979 : END IF
3980 410 : IF (do_ener) THEN
3981 208 : DO icoef = 1, SIZE(qs_ot_env(ichannel)%ener_x)
3982 : qs_ot_env(ichannel)%ener_x(icoef) = SIN(REAL(icoef + &
3983 : 37*qs_ot_env(ichannel)%spin_index + &
3984 208 : 101*qs_ot_env(ichannel)%kpoint_index, KIND=dp))
3985 : END DO
3986 : END IF
3987 : END DO
3988 192 : END SUBROUTINE broyden_randomize_current
3989 :
3990 : ! **************************************************************************************************
3991 : !> \brief Set all current product coordinates to zero.
3992 : ! **************************************************************************************************
3993 192 : SUBROUTINE broyden_set_current_zero()
3994 : INTEGER :: ichannel
3995 :
3996 410 : DO ichannel = 1, nspin
3997 218 : CALL dbcsr_set(qs_ot_env(ichannel)%matrix_x, 0.0_dp)
3998 218 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
3999 52 : CALL dbcsr_set(qs_ot_env(ichannel)%matrix_x_im, 0.0_dp)
4000 : END IF
4001 218 : IF (do_rotation) THEN
4002 52 : CALL dbcsr_set(qs_ot_env(ichannel)%rot_mat_x, 0.0_dp)
4003 52 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
4004 52 : CALL dbcsr_set(qs_ot_env(ichannel)%rot_mat_x_im, 0.0_dp)
4005 : END IF
4006 : END IF
4007 566 : IF (do_ener) qs_ot_env(ichannel)%ener_x = 0.0_dp
4008 : END DO
4009 192 : END SUBROUTINE broyden_set_current_zero
4010 :
4011 : ! **************************************************************************************************
4012 : !> \brief Add one stored product vector to the current coordinate.
4013 : !> \param vector_kind position or residual history
4014 : !> \param history_index source history slot
4015 : !> \param alpha_current scale of the current coordinate
4016 : !> \param beta_history scale of the stored coordinate
4017 : ! **************************************************************************************************
4018 1970 : SUBROUTINE broyden_add_history(vector_kind, history_index, alpha_current, beta_history)
4019 : INTEGER, INTENT(IN) :: vector_kind, history_index
4020 : REAL(KIND=dp), INTENT(IN) :: alpha_current, beta_history
4021 :
4022 : INTEGER :: ichannel
4023 : TYPE(dbcsr_type), POINTER :: source
4024 :
4025 4222 : DO ichannel = 1, nspin
4026 2252 : CALL broyden_matrix_pointer(vector_kind, history_index, ichannel, .FALSE., .FALSE., source)
4027 2252 : CALL dbcsr_add(qs_ot_env(ichannel)%matrix_x, source, alpha_current, beta_history)
4028 2252 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
4029 564 : CALL broyden_matrix_pointer(vector_kind, history_index, ichannel, .TRUE., .FALSE., source)
4030 564 : CALL dbcsr_add(qs_ot_env(ichannel)%matrix_x_im, source, alpha_current, beta_history)
4031 : END IF
4032 2252 : IF (do_rotation) THEN
4033 564 : CALL broyden_matrix_pointer(vector_kind, history_index, ichannel, .FALSE., .TRUE., source)
4034 564 : CALL dbcsr_add(qs_ot_env(ichannel)%rot_mat_x, source, alpha_current, beta_history)
4035 564 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
4036 564 : CALL broyden_matrix_pointer(vector_kind, history_index, ichannel, .TRUE., .TRUE., source)
4037 564 : CALL dbcsr_add(qs_ot_env(ichannel)%rot_mat_x_im, source, alpha_current, beta_history)
4038 : END IF
4039 : END IF
4040 4222 : IF (do_ener) THEN
4041 872 : SELECT CASE (vector_kind)
4042 : CASE (broyden_position)
4043 : qs_ot_env(ichannel)%ener_x = alpha_current*qs_ot_env(ichannel)%ener_x + &
4044 1232 : beta_history*qs_ot_env(ichannel)%ener_h_x(history_index, :)
4045 : CASE (broyden_residual)
4046 : qs_ot_env(ichannel)%ener_x = alpha_current*qs_ot_env(ichannel)%ener_x + &
4047 1024 : beta_history*qs_ot_env(ichannel)%ener_h_e(history_index, :)
4048 : CASE DEFAULT
4049 564 : CPABORT("Invalid Broyden history vector kind")
4050 : END SELECT
4051 : END IF
4052 : END DO
4053 1970 : END SUBROUTINE broyden_add_history
4054 :
4055 : ! **************************************************************************************************
4056 : !> \brief Product-space scalar product including complex, rotation, and energy coordinates.
4057 : !> \param kind_a first vector kind
4058 : !> \param index_a first history index, ignored for current/gradient vectors
4059 : !> \param kind_b second vector kind
4060 : !> \param index_b second history index, ignored for current/gradient vectors
4061 : !> \param value globally reduced scalar product
4062 : ! **************************************************************************************************
4063 15604 : SUBROUTINE broyden_product_dot(kind_a, index_a, kind_b, index_b, value)
4064 : INTEGER, INTENT(IN) :: kind_a, index_a, kind_b, index_b
4065 : REAL(KIND=dp), INTENT(OUT) :: value
4066 :
4067 : INTEGER :: ichannel
4068 : REAL(KIND=dp) :: dot_value
4069 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_a, energy_b
4070 : TYPE(dbcsr_type), POINTER :: matrix_a, matrix_b
4071 :
4072 15604 : value = 0.0_dp
4073 33124 : DO ichannel = 1, nspin
4074 17520 : CALL broyden_matrix_pointer(kind_a, index_a, ichannel, .FALSE., .FALSE., matrix_a)
4075 17520 : CALL broyden_matrix_pointer(kind_b, index_b, ichannel, .FALSE., .FALSE., matrix_b)
4076 17520 : CALL dbcsr_dot(matrix_a, matrix_b, dot_value)
4077 17520 : value = value + dot_value
4078 17520 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
4079 3832 : CALL broyden_matrix_pointer(kind_a, index_a, ichannel, .TRUE., .FALSE., matrix_a)
4080 3832 : CALL broyden_matrix_pointer(kind_b, index_b, ichannel, .TRUE., .FALSE., matrix_b)
4081 3832 : CALL dbcsr_dot(matrix_a, matrix_b, dot_value)
4082 3832 : value = value + dot_value
4083 : END IF
4084 17520 : IF (do_rotation) THEN
4085 3832 : CALL broyden_matrix_pointer(kind_a, index_a, ichannel, .FALSE., .TRUE., matrix_a)
4086 3832 : CALL broyden_matrix_pointer(kind_b, index_b, ichannel, .FALSE., .TRUE., matrix_b)
4087 3832 : CALL dbcsr_dot(matrix_a, matrix_b, dot_value)
4088 3832 : value = value + 0.5_dp*dot_value
4089 3832 : IF (qs_ot_env(ichannel)%has_complex_kpoint_state) THEN
4090 3832 : CALL broyden_matrix_pointer(kind_a, index_a, ichannel, .TRUE., .TRUE., matrix_a)
4091 3832 : CALL broyden_matrix_pointer(kind_b, index_b, ichannel, .TRUE., .TRUE., matrix_b)
4092 3832 : CALL dbcsr_dot(matrix_a, matrix_b, dot_value)
4093 3832 : value = value + 0.5_dp*dot_value
4094 : END IF
4095 : END IF
4096 50644 : IF (do_ener) THEN
4097 3832 : CALL broyden_energy_pointer(kind_a, index_a, ichannel, energy_a)
4098 3832 : CALL broyden_energy_pointer(kind_b, index_b, ichannel, energy_b)
4099 15328 : value = value + DOT_PRODUCT(energy_a, energy_b)
4100 : END IF
4101 : END DO
4102 15604 : CALL ot_mini_sum_kpoint_scalar(para_env_inter_kp, value)
4103 15604 : END SUBROUTINE broyden_product_dot
4104 :
4105 : ! **************************************************************************************************
4106 : !> \brief Select one orbital or rotation matrix from the product-vector storage.
4107 : !> \param vector_kind vector kind
4108 : !> \param history_index history slot
4109 : !> \param channel local spin/k-point channel
4110 : !> \param imaginary select imaginary component
4111 : !> \param rotation select rotation instead of orbital component
4112 : !> \param matrix selected matrix
4113 : ! **************************************************************************************************
4114 61976 : SUBROUTINE broyden_matrix_pointer(vector_kind, history_index, channel, imaginary, rotation, matrix)
4115 : INTEGER, INTENT(IN) :: vector_kind, history_index, channel
4116 : LOGICAL, INTENT(IN) :: imaginary, rotation
4117 : TYPE(dbcsr_type), POINTER :: matrix
4118 :
4119 61976 : IF (imaginary) THEN
4120 16456 : CPASSERT(qs_ot_env(channel)%has_complex_kpoint_state)
4121 : END IF
4122 91610 : SELECT CASE (vector_kind)
4123 : CASE (broyden_position)
4124 29634 : CPASSERT(history_index > 0)
4125 29634 : IF (rotation) THEN
4126 7656 : IF (imaginary) THEN
4127 3828 : matrix => qs_ot_env(channel)%rot_mat_h_x_im(history_index)%matrix
4128 : ELSE
4129 3828 : matrix => qs_ot_env(channel)%rot_mat_h_x(history_index)%matrix
4130 : END IF
4131 : ELSE
4132 21978 : IF (imaginary) THEN
4133 3828 : matrix => qs_ot_env(channel)%matrix_h_x_im(history_index)%matrix
4134 : ELSE
4135 18150 : matrix => qs_ot_env(channel)%matrix_h_x(history_index)%matrix
4136 : END IF
4137 : END IF
4138 : CASE (broyden_residual)
4139 29858 : CPASSERT(history_index > 0)
4140 29858 : IF (rotation) THEN
4141 7864 : IF (imaginary) THEN
4142 3932 : matrix => qs_ot_env(channel)%rot_mat_h_e_im(history_index)%matrix
4143 : ELSE
4144 3932 : matrix => qs_ot_env(channel)%rot_mat_h_e(history_index)%matrix
4145 : END IF
4146 : ELSE
4147 21994 : IF (imaginary) THEN
4148 3932 : matrix => qs_ot_env(channel)%matrix_h_e_im(history_index)%matrix
4149 : ELSE
4150 18062 : matrix => qs_ot_env(channel)%matrix_h_e(history_index)%matrix
4151 : END IF
4152 : END IF
4153 : CASE (broyden_random)
4154 1928 : IF (rotation) THEN
4155 728 : IF (imaginary) THEN
4156 364 : matrix => qs_ot_env(channel)%rot_mat_x_im
4157 : ELSE
4158 364 : matrix => qs_ot_env(channel)%rot_mat_x
4159 : END IF
4160 : ELSE
4161 1200 : IF (imaginary) THEN
4162 364 : matrix => qs_ot_env(channel)%matrix_x_im
4163 : ELSE
4164 836 : matrix => qs_ot_env(channel)%matrix_x
4165 : END IF
4166 : END IF
4167 : CASE (broyden_gradient)
4168 556 : IF (rotation) THEN
4169 208 : IF (imaginary) THEN
4170 104 : matrix => qs_ot_env(channel)%rot_mat_gx_im
4171 : ELSE
4172 104 : matrix => qs_ot_env(channel)%rot_mat_gx
4173 : END IF
4174 : ELSE
4175 348 : IF (imaginary) THEN
4176 104 : matrix => qs_ot_env(channel)%matrix_gx_im
4177 : ELSE
4178 244 : matrix => qs_ot_env(channel)%matrix_gx
4179 : END IF
4180 : END IF
4181 : CASE DEFAULT
4182 61976 : CPABORT("Invalid Broyden product vector kind")
4183 : END SELECT
4184 61976 : END SUBROUTINE broyden_matrix_pointer
4185 :
4186 : ! **************************************************************************************************
4187 : !> \brief Select one energy-variable vector from the product-vector storage.
4188 : !> \param vector_kind vector kind
4189 : !> \param history_index history slot
4190 : !> \param channel local spin/k-point channel
4191 : !> \param energy selected energy vector
4192 : ! **************************************************************************************************
4193 7664 : SUBROUTINE broyden_energy_pointer(vector_kind, history_index, channel, energy)
4194 : INTEGER, INTENT(IN) :: vector_kind, history_index, channel
4195 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy
4196 :
4197 11184 : SELECT CASE (vector_kind)
4198 : CASE (broyden_position)
4199 3520 : CPASSERT(history_index > 0)
4200 3520 : energy => qs_ot_env(channel)%ener_h_x(history_index, :)
4201 : CASE (broyden_residual)
4202 3676 : CPASSERT(history_index > 0)
4203 3676 : energy => qs_ot_env(channel)%ener_h_e(history_index, :)
4204 : CASE (broyden_random)
4205 364 : energy => qs_ot_env(channel)%ener_x
4206 : CASE (broyden_gradient)
4207 104 : energy => qs_ot_env(channel)%ener_gx
4208 : CASE DEFAULT
4209 7664 : CPABORT("Invalid Broyden energy vector kind")
4210 : END SELECT
4211 7664 : END SUBROUTINE broyden_energy_pointer
4212 :
4213 : END SUBROUTINE ot_broyden_step
4214 :
4215 : ! **************************************************************************************************
4216 : !> \brief ...
4217 : !> \param G ...
4218 : !> \param S ...
4219 : !> \param n ...
4220 : !> \return ...
4221 : ! **************************************************************************************************
4222 166 : FUNCTION new_sigma(G, S, n) RESULT(sigma)
4223 : !
4224 : ! Calculate new sigma from eigenvalues of full size G by Arnoldi.
4225 : !
4226 : ! **************************************************************************************************
4227 :
4228 : REAL(KIND=dp), DIMENSION(:, :) :: G, S
4229 : INTEGER :: n
4230 : REAL(KIND=dp) :: sigma
4231 :
4232 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigv
4233 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H
4234 :
4235 664 : ALLOCATE (H(n, n))
4236 166 : CALL hess_G(G, S, H, n)
4237 498 : ALLOCATE (eigv(n))
4238 166 : CALL diamat_all(H(1:n, 1:n), eigv)
4239 :
4240 : SELECT CASE (1)
4241 : CASE (1)
4242 : ! This estimator seems to work well. No theory.
4243 2064 : sigma = SUM(ABS(eigv**2))/SUM(ABS(eigv))
4244 : CASE (2)
4245 : ! Estimator based on Frobenius norm minimizer
4246 : sigma = SUM(ABS(eigv))/MAX(1, SIZE(eigv))
4247 : CASE (3)
4248 : ! Estimator based on induced 2-norm
4249 : sigma = (MAXVAL(ABS(eigv)) + MINVAL(ABS(eigv)))*0.5_dp
4250 : END SELECT
4251 :
4252 166 : DEALLOCATE (H, eigv)
4253 166 : END FUNCTION new_sigma
4254 :
4255 : ! **************************************************************************************************
4256 : !> \brief ...
4257 : !> \param G ...
4258 : !> \param S ...
4259 : !> \param H ...
4260 : !> \param n ...
4261 : ! **************************************************************************************************
4262 166 : SUBROUTINE hess_G(G, S, H, n)
4263 : !
4264 : ! Make a hessenberg out of G into H. Cf Arnoldi.
4265 : ! Inner product is weighted by S.
4266 : ! Possible lucky breakdown at n.
4267 : !
4268 : ! **************************************************************************************************
4269 : REAL(KIND=dp), DIMENSION(:, :) :: G, S, H
4270 : INTEGER :: n
4271 :
4272 : INTEGER :: i, j, k
4273 : REAL(KIND=dp) :: tmp
4274 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: v
4275 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Q
4276 :
4277 166 : i = SIZE(G, 1)
4278 166 : k = SIZE(H, 1)
4279 664 : ALLOCATE (Q(i, k))
4280 498 : ALLOCATE (v(i))
4281 6458 : H = 0.0_dp
4282 166 : Q = 0.0_dp
4283 :
4284 2084 : Q(:, 1) = 1.0_dp
4285 29668 : tmp = SQRT(DOT_PRODUCT(Q(:, 1), MATMUL(S, Q(:, 1))))
4286 12750 : Q(:, :) = Q(:, :)/tmp
4287 :
4288 1028 : DO i = 1, k
4289 174698 : v(:) = MATMUL(G, Q(:, i))
4290 3954 : DO j = 1, i
4291 715800 : H(j, i) = DOT_PRODUCT(Q(:, j), MATMUL(S, v))
4292 47142 : v(:) = v - H(j, i)*Q(:, j)
4293 : END DO
4294 1028 : IF (i < k) THEN
4295 159352 : tmp = DOT_PRODUCT(v, MATMUL(S, v))
4296 704 : IF (tmp <= 0.0_dp) THEN
4297 0 : n = i
4298 0 : EXIT
4299 : END IF
4300 704 : tmp = SQRT(tmp)
4301 : ! Lucky breakdown
4302 704 : IF (ABS(tmp) < 1e-9_dp) THEN
4303 4 : n = i
4304 4 : EXIT
4305 : END IF
4306 700 : H(i + 1, i) = tmp
4307 10340 : Q(:, i + 1) = v/H(i + 1, i)
4308 : END IF
4309 : END DO
4310 :
4311 166 : DEALLOCATE (Q, v)
4312 166 : END SUBROUTINE hess_G
4313 :
4314 6088 : END MODULE qs_ot_minimizer
|
