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