(git:b279b6b)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief orbital transformations
10 !> \par History
11 !> None
12 !> \author Joost VandeVondele (09.2002)
13 ! **************************************************************************************************
14 MODULE qs_ot_minimizer
15 
16  USE cp_log_handling, ONLY: cp_get_default_logger,&
17  cp_logger_get_default_unit_nr,&
18  cp_logger_type
19  USE dbcsr_api, ONLY: dbcsr_add,&
20  dbcsr_copy,&
21  dbcsr_dot,&
22  dbcsr_get_info,&
23  dbcsr_init_random,&
24  dbcsr_p_type,&
25  dbcsr_scale,&
26  dbcsr_set
27  USE kinds, ONLY: dp,&
28  int_8
29  USE mathlib, ONLY: diamat_all
30  USE preconditioner, ONLY: apply_preconditioner
31  USE qs_ot, ONLY: qs_ot_get_derivative,&
32  qs_ot_get_derivative_ref
33  USE qs_ot_types, ONLY: qs_ot_type
34 #include "./base/base_uses.f90"
35 
36  IMPLICIT NONE
37 
38  PRIVATE
39 
40  PUBLIC :: ot_mini
41 
42  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_minimizer'
43 
44 CONTAINS
45 !
46 ! the minimizer interface
47 ! should present all possible modes of minimization
48 ! these include CG SD DIIS
49 !
50 !
51 ! IN the case of nspin != 1 we have a gradient that is distributed over different qs_ot_env.
52 ! still things remain basically the same, since there are no constraints between the different qs_ot_env
53 ! we only should take care that the various scalar products are taken over the full vectors.
54 ! all the information needed and collected can be stored in the fist qs_ot_env only
55 ! (indicating that the data type for the gradient/position and minization should be separated)
56 !
57 ! **************************************************************************************************
58 !> \brief ...
59 !> \param qs_ot_env ...
60 !> \param matrix_hc ...
61 ! **************************************************************************************************
62  SUBROUTINE ot_mini(qs_ot_env, matrix_hc)
63  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
64  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc
65 
66  CHARACTER(len=*), PARAMETER :: routinen = 'ot_mini'
67 
68  INTEGER :: handle, ispin, nspin
69  LOGICAL :: do_ener, do_ks
70  REAL(kind=dp) :: tmp
71 
72  CALL timeset(routinen, handle)
73 
74  nspin = SIZE(qs_ot_env)
75 
76  do_ks = qs_ot_env(1)%settings%ks
77  do_ener = qs_ot_env(1)%settings%do_ener
78 
79  qs_ot_env(1)%OT_METHOD_FULL = ""
80 
81  ! compute the gradient for the variables x
82  IF (.NOT. qs_ot_env(1)%energy_only) THEN
83  qs_ot_env(1)%gradient = 0.0_dp
84  DO ispin = 1, nspin
85  IF (do_ks) THEN
86  SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
87  CASE ("TOD")
88  CALL qs_ot_get_derivative(matrix_hc(ispin)%matrix, qs_ot_env(ispin)%matrix_x, &
89  qs_ot_env(ispin)%matrix_sx, &
90  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
91  CASE ("REF")
92  CALL qs_ot_get_derivative_ref(matrix_hc(ispin)%matrix, &
93  qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, &
94  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
95  CASE DEFAULT
96  cpabort("ALGORITHM NYI")
97  END SELECT
98  END IF
99  ! and also the gradient along the direction
100  IF (qs_ot_env(1)%use_dx) THEN
101  IF (do_ks) THEN
102  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
103  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
104  IF (qs_ot_env(1)%settings%do_rotation) THEN
105  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
106  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + 0.5_dp*tmp
107  END IF
108  END IF
109  IF (do_ener) THEN
110  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
111  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
112  END IF
113  ELSE
114  IF (do_ks) THEN
115  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
116  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
117  IF (qs_ot_env(1)%settings%do_rotation) THEN
118  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
119  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - 0.5_dp*tmp
120  END IF
121  END IF
122  IF (do_ener) THEN
123  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
124  qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
125  END IF
126  END IF
127  END DO
128  END IF
129 
130  SELECT CASE (qs_ot_env(1)%settings%OT_METHOD)
131  CASE ("CG")
132  IF (current_point_is_fine(qs_ot_env)) THEN
133  qs_ot_env(1)%OT_METHOD_FULL = "OT CG"
134  CALL ot_new_cg_direction(qs_ot_env)
135  qs_ot_env(1)%line_search_count = 0
136  ELSE
137  qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
138  END IF
139  CALL do_line_search(qs_ot_env)
140  CASE ("SD")
141  IF (current_point_is_fine(qs_ot_env)) THEN
142  qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
143  CALL ot_new_sd_direction(qs_ot_env)
144  qs_ot_env(1)%line_search_count = 0
145  ELSE
146  qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
147  END IF
148  CALL do_line_search(qs_ot_env)
149  CASE ("DIIS")
150  qs_ot_env(1)%OT_METHOD_FULL = "OT DIIS"
151  CALL ot_diis_step(qs_ot_env)
152  CASE ("BROY")
153  qs_ot_env(1)%OT_METHOD_FULL = "OT BROY"
154  CALL ot_broyden_step(qs_ot_env)
155  CASE DEFAULT
156  cpabort("OT_METHOD NYI")
157  END SELECT
158 
159  CALL timestop(handle)
160 
161  END SUBROUTINE ot_mini
162 
163 !
164 ! checks if the current point is a good point for finding a new direction
165 ! or if we should improve the line_search, if it is used
166 !
167 ! **************************************************************************************************
168 !> \brief ...
169 !> \param qs_ot_env ...
170 !> \return ...
171 ! **************************************************************************************************
172  FUNCTION current_point_is_fine(qs_ot_env) RESULT(res)
173  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
174  LOGICAL :: res
175 
176  res = .false.
177 
178  ! only if we have a gradient it can be fine
179  IF (.NOT. qs_ot_env(1)%energy_only) THEN
180 
181  ! we have not yet started with the line search
182  IF (qs_ot_env(1)%line_search_count .EQ. 0) THEN
183  res = .true.
184  RETURN
185  END IF
186 
187  IF (qs_ot_env(1)%line_search_might_be_done) THEN
188  ! here we put the more complicated logic later
189  res = .true.
190  RETURN
191  END IF
192 
193  END IF
194 
195  END FUNCTION current_point_is_fine
196 
197 !
198 ! performs various kinds of line searches
199 !
200 ! **************************************************************************************************
201 !> \brief ...
202 !> \param qs_ot_env ...
203 ! **************************************************************************************************
204  SUBROUTINE do_line_search(qs_ot_env)
205  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
206 
207  SELECT CASE (qs_ot_env(1)%settings%line_search_method)
208  CASE ("GOLD")
209  CALL do_line_search_gold(qs_ot_env)
210  CASE ("3PNT")
211  CALL do_line_search_3pnt(qs_ot_env)
212  CASE ("2PNT")
213  CALL do_line_search_2pnt(qs_ot_env)
214  CASE ("NONE")
215  CALL do_line_search_none(qs_ot_env)
216  CASE DEFAULT
217  cpabort("NYI")
218  END SELECT
219  END SUBROUTINE do_line_search
220 
221 ! **************************************************************************************************
222 !> \brief moves x adding the right amount (ds) of the gradient or search direction
223 !> \param ds ...
224 !> \param qs_ot_env ...
225 !> \par History
226 !> 08.2004 created [ Joost VandeVondele ] copied here from a larger number of subroutines
227 ! **************************************************************************************************
228  SUBROUTINE take_step(ds, qs_ot_env)
229  REAL(kind=dp) :: ds
230  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
231 
232  CHARACTER(len=*), PARAMETER :: routinen = 'take_step'
233 
234  INTEGER :: handle, ispin, nspin
235  LOGICAL :: do_ener, do_ks
236 
237  CALL timeset(routinen, handle)
238 
239  nspin = SIZE(qs_ot_env)
240 
241  do_ks = qs_ot_env(1)%settings%ks
242  do_ener = qs_ot_env(1)%settings%do_ener
243 
244  ! now update x to take into account this new step
245  ! either dx or -gx is the direction to use
246  IF (qs_ot_env(1)%use_dx) THEN
247  IF (do_ks) THEN
248  DO ispin = 1, nspin
249  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_dx, &
250  alpha_scalar=1.0_dp, beta_scalar=ds)
251  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
252  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, qs_ot_env(ispin)%rot_mat_dx, &
253  alpha_scalar=1.0_dp, beta_scalar=ds)
254  END IF
255  END DO
256  END IF
257  IF (do_ener) THEN
258  DO ispin = 1, nspin
259  qs_ot_env(ispin)%ener_x = qs_ot_env(ispin)%ener_x + ds*qs_ot_env(ispin)%ener_dx
260  END DO
261  END IF
262  ELSE
263  IF (do_ks) THEN
264  DO ispin = 1, nspin
265  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_gx, &
266  alpha_scalar=1.0_dp, beta_scalar=-ds)
267  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
268  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, qs_ot_env(ispin)%rot_mat_gx, &
269  alpha_scalar=1.0_dp, beta_scalar=-ds)
270  END IF
271  END DO
272  END IF
273  IF (do_ener) THEN
274  DO ispin = 1, nspin
275  qs_ot_env(ispin)%ener_x = qs_ot_env(ispin)%ener_x - ds*qs_ot_env(ispin)%ener_gx
276  END DO
277  END IF
278  END IF
279  CALL timestop(handle)
280  END SUBROUTINE take_step
281 
282 ! implements a golden ratio search as a robust way of minimizing
283 ! **************************************************************************************************
284 !> \brief ...
285 !> \param qs_ot_env ...
286 ! **************************************************************************************************
287  SUBROUTINE do_line_search_gold(qs_ot_env)
288 
289  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
290 
291  CHARACTER(len=*), PARAMETER :: routinen = 'do_line_search_gold'
292  REAL(kind=dp), PARAMETER :: gold_sec = 0.3819_dp
293 
294  INTEGER :: count, handle
295  REAL(kind=dp) :: ds
296 
297  CALL timeset(routinen, handle)
298 
299  qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
300  count = qs_ot_env(1)%line_search_count
301  qs_ot_env(1)%line_search_might_be_done = .false.
302  qs_ot_env(1)%energy_only = .true.
303 
304  IF (count + 1 .GT. SIZE(qs_ot_env(1)%OT_pos)) THEN
305  ! should not happen, we pass with a warning first
306  ! you can increase the size of OT_pos and the like in qs_ot_env
307  cpabort("MAX ITER EXCEEDED : FATAL")
308  END IF
309 
310  IF (qs_ot_env(1)%line_search_count .EQ. 1) THEN
311  qs_ot_env(1)%line_search_left = 1
312  qs_ot_env(1)%line_search_right = 0
313  qs_ot_env(1)%line_search_mid = 1
314  qs_ot_env(1)%ot_pos(1) = 0.0_dp
315  qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
316  qs_ot_env(1)%ot_pos(2) = qs_ot_env(1)%ds_min/gold_sec
317  ELSE
318  qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
319  ! it's essentially a book keeping game.
320  ! keep left on the left, keep (bring) right on the right
321  ! and mid in between these two
322  IF (qs_ot_env(1)%line_search_right .EQ. 0) THEN ! we do not yet have the right bracket
323  IF (qs_ot_env(1)%ot_energy(count - 1) .LT. qs_ot_env(1)%ot_energy(count)) THEN
324  qs_ot_env(1)%line_search_right = count
325  qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) + &
326  (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) - &
327  qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))*gold_sec
328  ELSE
329  qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
330  qs_ot_env(1)%line_search_mid = count
331  qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ot_pos(count)/gold_sec ! expand
332  END IF
333  ELSE
334  ! first determine where we are and construct the new triplet
335  IF (qs_ot_env(1)%ot_pos(count) .LT. qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
336  IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
337  qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
338  qs_ot_env(1)%line_search_mid = count
339  ELSE
340  qs_ot_env(1)%line_search_left = count
341  END IF
342  ELSE
343  IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
344  qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
345  qs_ot_env(1)%line_search_mid = count
346  ELSE
347  qs_ot_env(1)%line_search_right = count
348  END IF
349  END IF
350  ! now find the new point in the largest section
351  IF ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
352  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) .GT. &
353  (qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
354  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left))) THEN
355  qs_ot_env(1)%ot_pos(count + 1) = &
356  qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) + &
357  gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
358  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))
359  ELSE
360  qs_ot_env(1)%ot_pos(count + 1) = &
361  qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left) + &
362  gold_sec*(qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
363  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left))
364  END IF
365  ! check for termination
366  IF (((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_right) &
367  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) .LT. &
368  qs_ot_env(1)%ds_min*qs_ot_env(1)%settings%gold_target) .AND. &
369  ((qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid) &
370  - qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_left)) .LT. &
371  qs_ot_env(1)%ds_min*qs_ot_env(1)%settings%gold_target)) THEN
372  qs_ot_env(1)%energy_only = .false.
373  qs_ot_env(1)%line_search_might_be_done = .true.
374  END IF
375  END IF
376  END IF
377  ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
378  qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
379 
380  CALL take_step(ds, qs_ot_env)
381 
382  CALL timestop(handle)
383 
384  END SUBROUTINE do_line_search_gold
385 
386 ! **************************************************************************************************
387 !> \brief ...
388 !> \param qs_ot_env ...
389 ! **************************************************************************************************
390  SUBROUTINE do_line_search_3pnt(qs_ot_env)
391 
392  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
393 
394  CHARACTER(len=*), PARAMETER :: routinen = 'do_line_search_3pnt'
395 
396  INTEGER :: count, handle
397  REAL(kind=dp) :: denom, ds, fa, fb, fc, nom, pos, val, &
398  xa, xb, xc
399 
400  CALL timeset(routinen, handle)
401 
402  qs_ot_env(1)%line_search_might_be_done = .false.
403  qs_ot_env(1)%energy_only = .true.
404 
405  ! a three point interpolation based on the energy
406  qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
407  count = qs_ot_env(1)%line_search_count
408  qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
409  SELECT CASE (count)
410  CASE (1)
411  qs_ot_env(1)%ot_pos(count) = 0.0_dp
412  qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*0.8_dp
413  CASE (2)
414  IF (qs_ot_env(1)%OT_energy(count) .GT. qs_ot_env(1)%OT_energy(count - 1)) THEN
415  qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*0.5_dp
416  ELSE
417  qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*1.4_dp
418  END IF
419  CASE (3)
420  xa = qs_ot_env(1)%OT_pos(1)
421  xb = qs_ot_env(1)%OT_pos(2)
422  xc = qs_ot_env(1)%OT_pos(3)
423  fa = qs_ot_env(1)%OT_energy(1)
424  fb = qs_ot_env(1)%OT_energy(2)
425  fc = qs_ot_env(1)%OT_energy(3)
426  nom = (xb - xa)**2*(fb - fc) - (xb - xc)**2*(fb - fa)
427  denom = (xb - xa)*(fb - fc) - (xb - xc)*(fb - fa)
428  IF (abs(denom) .LE. 1.0e-18_dp*max(abs(fb - fc), abs(fb - fa))) THEN
429  pos = xb
430  ELSE
431  pos = xb - 0.5_dp*nom/denom ! position of the stationary point
432  END IF
433  val = (pos - xa)*(pos - xb)*fc/((xc - xa)*(xc - xb)) + &
434  (pos - xb)*(pos - xc)*fa/((xa - xb)*(xa - xc)) + &
435  (pos - xc)*(pos - xa)*fb/((xb - xc)*(xb - xa))
436  IF (val .LT. fa .AND. val .LE. fb .AND. val .LE. fc) THEN ! OK, we go to a minimum
437  ! we take a guard against too large steps
438  qs_ot_env(1)%OT_pos(count + 1) = max(maxval(qs_ot_env(1)%OT_pos(1:3))*0.01_dp, &
439  min(pos, maxval(qs_ot_env(1)%OT_pos(1:3))*4.0_dp))
440  ELSE ! just take an extended step
441  qs_ot_env(1)%OT_pos(count + 1) = maxval(qs_ot_env(1)%OT_pos(1:3))*2.0_dp
442  END IF
443  qs_ot_env(1)%energy_only = .false.
444  qs_ot_env(1)%line_search_might_be_done = .true.
445  CASE DEFAULT
446  cpabort("NYI")
447  END SELECT
448  ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
449  qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
450 
451  CALL take_step(ds, qs_ot_env)
452 
453  CALL timestop(handle)
454 
455  END SUBROUTINE do_line_search_3pnt
456 
457 ! **************************************************************************************************
458 !> \brief ...
459 !> \param qs_ot_env ...
460 ! **************************************************************************************************
461  SUBROUTINE do_line_search_2pnt(qs_ot_env)
462 
463  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
464 
465  CHARACTER(len=*), PARAMETER :: routinen = 'do_line_search_2pnt'
466 
467  INTEGER :: count, handle
468  REAL(kind=dp) :: a, b, c, ds, pos, val, x0, x1
469 
470  CALL timeset(routinen, handle)
471 
472  qs_ot_env(1)%line_search_might_be_done = .false.
473  qs_ot_env(1)%energy_only = .true.
474 
475  ! a three point interpolation based on the energy
476  qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
477  count = qs_ot_env(1)%line_search_count
478  qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
479  SELECT CASE (count)
480  CASE (1)
481  qs_ot_env(1)%ot_pos(count) = 0.0_dp
482  qs_ot_env(1)%ot_grad(count) = qs_ot_env(1)%gradient
483  qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*1.0_dp
484  CASE (2)
485  x0 = 0.0_dp
486  c = qs_ot_env(1)%ot_energy(1)
487  b = qs_ot_env(1)%ot_grad(1)
488  x1 = qs_ot_env(1)%ot_pos(2)
489  a = (qs_ot_env(1)%ot_energy(2) - b*x1 - c)/(x1**2)
490  IF (a .LE. 0.0_dp) a = 1.0e-15_dp
491  pos = -b/(2.0_dp*a)
492  val = a*pos**2 + b*pos + c
493  qs_ot_env(1)%energy_only = .false.
494  qs_ot_env(1)%line_search_might_be_done = .true.
495  IF (val .LT. qs_ot_env(1)%ot_energy(1) .AND. val .LE. qs_ot_env(1)%ot_energy(2)) THEN
496  ! we go to a minimum, but ...
497  ! we take a guard against too large steps
498  qs_ot_env(1)%OT_pos(count + 1) = max(maxval(qs_ot_env(1)%OT_pos(1:2))*0.01_dp, &
499  min(pos, maxval(qs_ot_env(1)%OT_pos(1:2))*4.0_dp))
500  ELSE ! just take an extended step
501  qs_ot_env(1)%OT_pos(count + 1) = maxval(qs_ot_env(1)%OT_pos(1:2))*2.0_dp
502  END IF
503  CASE DEFAULT
504  cpabort("NYI")
505  END SELECT
506  ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
507  qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
508 
509  CALL take_step(ds, qs_ot_env)
510 
511  CALL timestop(handle)
512 
513  END SUBROUTINE do_line_search_2pnt
514 
515 ! **************************************************************************************************
516 !> \brief ...
517 !> \param qs_ot_env ...
518 ! **************************************************************************************************
519  SUBROUTINE do_line_search_none(qs_ot_env)
520  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
521 
522  CALL take_step(qs_ot_env(1)%ds_min, qs_ot_env)
523 
524  END SUBROUTINE do_line_search_none
525 
526 !
527 ! creates a new SD direction, using the preconditioner if associated
528 ! also updates the gradient for line search
529 !
530 
531 ! **************************************************************************************************
532 !> \brief ...
533 !> \param qs_ot_env ...
534 ! **************************************************************************************************
535  SUBROUTINE ot_new_sd_direction(qs_ot_env)
536  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
537 
538  CHARACTER(len=*), PARAMETER :: routinen = 'ot_new_sd_direction'
539 
540  INTEGER :: handle, ispin, itmp, k, n, nener, nspin
541  LOGICAL :: do_ener, do_ks
542  REAL(kind=dp) :: tmp
543  TYPE(cp_logger_type), POINTER :: logger
544 
545  CALL timeset(routinen, handle)
546 
547 !***SCP
548 
549  nspin = SIZE(qs_ot_env)
550  logger => cp_get_default_logger()
551  do_ks = qs_ot_env(1)%settings%ks
552  do_ener = qs_ot_env(1)%settings%do_ener
553 
554  IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
555  IF (.NOT. qs_ot_env(1)%use_dx) cpabort("use dx")
556  qs_ot_env(1)%gnorm = 0.0_dp
557  IF (do_ks) THEN
558  DO ispin = 1, nspin
559  CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
560  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx)
561  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
562  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
563  END DO
564  IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
565  logger => cp_get_default_logger()
566  WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
567  END IF
568  DO ispin = 1, nspin
569  CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
570  END DO
571  IF (qs_ot_env(1)%settings%do_rotation) THEN
572  DO ispin = 1, nspin
573  ! right now no preconditioner yet
574  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
575  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
576  ! added 0.5, because we have (antisymmetry) only half the number of variables
577  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
578  END DO
579  DO ispin = 1, nspin
580  CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
581  END DO
582  END IF
583  END IF
584  IF (do_ener) THEN
585  DO ispin = 1, nspin
586  qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_gx
587  tmp = dot_product(qs_ot_env(ispin)%ener_dx, qs_ot_env(ispin)%ener_gx)
588  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
589  qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_dx
590  END DO
591  END IF
592  ELSE
593  qs_ot_env(1)%gnorm = 0.0_dp
594  IF (do_ks) THEN
595  DO ispin = 1, nspin
596  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
597  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
598  END DO
599  IF (qs_ot_env(1)%settings%do_rotation) THEN
600  DO ispin = 1, nspin
601  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
602  ! added 0.5, because we have (antisymmetry) only half the number of variables
603  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
604  END DO
605  END IF
606  END IF
607  IF (do_ener) THEN
608  DO ispin = 1, nspin
609  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
610  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
611  END DO
612  END IF
613  END IF
614 
615  k = 0
616  n = 0
617  nener = 0
618  IF (do_ks) THEN
619  CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
620  DO ispin = 1, nspin
621  CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
622  k = k + itmp
623  END DO
624  END IF
625  IF (do_ener) THEN
626  DO ispin = 1, nspin
627  nener = nener + SIZE(qs_ot_env(ispin)%ener_x)
628  END DO
629  END IF
630  ! Handling the case of no free variables to optimize
631  IF (int(n, kind=int_8)*int(k, kind=int_8) + nener /= 0) THEN
632  qs_ot_env(1)%delta = sqrt(abs(qs_ot_env(1)%gnorm)/(int(n, kind=int_8)*int(k, kind=int_8) + nener))
633  qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
634  ELSE
635  qs_ot_env(1)%delta = 0.0_dp
636  qs_ot_env(1)%gradient = 0.0_dp
637  END IF
638 
639  CALL timestop(handle)
640 
641  END SUBROUTINE ot_new_sd_direction
642 
643 !
644 ! creates a new CG direction. Implements Polak-Ribierre variant
645 ! using the preconditioner if associated
646 ! also updates the gradient for line search
647 !
648 ! **************************************************************************************************
649 !> \brief ...
650 !> \param qs_ot_env ...
651 ! **************************************************************************************************
652  SUBROUTINE ot_new_cg_direction(qs_ot_env)
653  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
654 
655  CHARACTER(len=*), PARAMETER :: routinen = 'ot_new_cg_direction'
656 
657  INTEGER :: handle, ispin, itmp, k, n, nener, nspin
658  LOGICAL :: do_ener, do_ks
659  REAL(kind=dp) :: beta_pr, gnorm_cross, test_down, tmp
660  TYPE(cp_logger_type), POINTER :: logger
661 
662  CALL timeset(routinen, handle)
663 
664  nspin = SIZE(qs_ot_env)
665  logger => cp_get_default_logger()
666 
667  do_ks = qs_ot_env(1)%settings%ks
668  do_ener = qs_ot_env(1)%settings%do_ener
669 
670  gnorm_cross = 0.0_dp
671  IF (do_ks) THEN
672  DO ispin = 1, nspin
673  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
674  gnorm_cross = gnorm_cross + tmp
675  END DO
676  IF (qs_ot_env(1)%settings%do_rotation) THEN
677  DO ispin = 1, nspin
678  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
679  ! added 0.5, because we have (antisymmetry) only half the number of variables
680  gnorm_cross = gnorm_cross + 0.5_dp*tmp
681  END DO
682  END IF
683  END IF
684  IF (do_ener) THEN
685  DO ispin = 1, nspin
686  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
687  gnorm_cross = gnorm_cross + tmp
688  END DO
689  END IF
690 
691  IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
692 
693  DO ispin = 1, nspin
694  CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
695  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
696  END DO
697  qs_ot_env(1)%gnorm = 0.0_dp
698  IF (do_ks) THEN
699  DO ispin = 1, nspin
700  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
701  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
702  END DO
703  IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
704  WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
705  END IF
706  DO ispin = 1, nspin
707  CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
708  END DO
709  IF (qs_ot_env(1)%settings%do_rotation) THEN
710  DO ispin = 1, nspin
711  ! right now no preconditioner yet
712  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
713  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old, tmp)
714  ! added 0.5, because we have (antisymmetry) only half the number of variables
715  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
716  END DO
717  DO ispin = 1, nspin
718  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx_old)
719  END DO
720  END IF
721  END IF
722  IF (do_ener) THEN
723  DO ispin = 1, nspin
724  qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
725  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
726  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
727  qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_gx_old
728  END DO
729  END IF
730  ELSE
731  IF (do_ks) THEN
732  qs_ot_env(1)%gnorm = 0.0_dp
733  DO ispin = 1, nspin
734  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
735  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
736  CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx)
737  END DO
738  IF (qs_ot_env(1)%settings%do_rotation) THEN
739  DO ispin = 1, nspin
740  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
741  ! added 0.5, because we have (antisymmetry) only half the number of variables
742  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
743  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
744  END DO
745  END IF
746  END IF
747  IF (do_ener) THEN
748  DO ispin = 1, nspin
749  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
750  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
751  qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
752  END DO
753  END IF
754  END IF
755 
756  k = 0
757  n = 0
758  nener = 0
759  IF (do_ks) THEN
760  CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
761  DO ispin = 1, nspin
762  CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
763  k = k + itmp
764  END DO
765  END IF
766  IF (do_ener) THEN
767  DO ispin = 1, nspin
768  nener = nener + SIZE(qs_ot_env(ispin)%ener_x)
769  END DO
770  END IF
771  ! Handling the case of no free variables to optimize
772  IF (int(n, kind=int_8)*int(k, kind=int_8) + nener /= 0) THEN
773  qs_ot_env(1)%delta = sqrt(abs(qs_ot_env(1)%gnorm)/(int(n, kind=int_8)*int(k, kind=int_8) + nener))
774  beta_pr = (qs_ot_env(1)%gnorm - gnorm_cross)/qs_ot_env(1)%gnorm_old
775  ELSE
776  qs_ot_env(1)%delta = 0.0_dp
777  beta_pr = 0.0_dp
778  END IF
779  beta_pr = max(beta_pr, 0.0_dp) ! reset to SD
780 
781  test_down = 0.0_dp
782  IF (do_ks) THEN
783  DO ispin = 1, nspin
784  CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
785  alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
786  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
787  test_down = test_down + tmp
788  IF (qs_ot_env(1)%settings%do_rotation) THEN
789  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx, &
790  alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
791  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
792  test_down = test_down + 0.5_dp*tmp
793  END IF
794  END DO
795  END IF
796  IF (do_ener) THEN
797  DO ispin = 1, nspin
798  qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - qs_ot_env(ispin)%ener_gx
799  tmp = dot_product(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
800  test_down = test_down + tmp
801  END DO
802  END IF
803 
804  IF (test_down .GE. 0.0_dp) THEN ! reset to SD
805  beta_pr = 0.0_dp
806  IF (do_ks) THEN
807  DO ispin = 1, nspin
808  CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
809  alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
810  IF (qs_ot_env(1)%settings%do_rotation) THEN
811  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
812  qs_ot_env(ispin)%rot_mat_gx, alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
813  END IF
814  END DO
815  END IF
816  IF (do_ener) THEN
817  DO ispin = 1, nspin
818  qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - qs_ot_env(ispin)%ener_gx
819  END DO
820  END IF
821  END IF
822  ! since we change the direction we have to adjust the gradient
823  qs_ot_env(1)%gradient = beta_pr*qs_ot_env(1)%gradient - qs_ot_env(1)%gnorm
824  qs_ot_env(1)%gnorm_old = qs_ot_env(1)%gnorm
825 
826  CALL timestop(handle)
827 
828  END SUBROUTINE ot_new_cg_direction
829 
830 ! **************************************************************************************************
831 !> \brief ...
832 !> \param qs_ot_env ...
833 ! **************************************************************************************************
834  SUBROUTINE ot_diis_step(qs_ot_env)
835  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
836 
837  CHARACTER(len=*), PARAMETER :: routinen = 'ot_diis_step'
838 
839  INTEGER :: diis_bound, diis_m, handle, i, info, &
840  ispin, itmp, j, k, n, nener, nspin
841  LOGICAL :: do_ener, do_ks, do_ot_sd
842  REAL(kind=dp) :: overlap, tmp, tr_xnew_gx, tr_xold_gx
843  TYPE(cp_logger_type), POINTER :: logger
844 
845  CALL timeset(routinen, handle)
846 
847  logger => cp_get_default_logger()
848 
849  do_ks = qs_ot_env(1)%settings%ks
850  do_ener = qs_ot_env(1)%settings%do_ener
851  nspin = SIZE(qs_ot_env)
852 
853  diis_m = qs_ot_env(1)%settings%diis_m
854 
855  IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
856  diis_bound = qs_ot_env(1)%diis_iter + 1
857  ELSE
858  diis_bound = diis_m
859  END IF
860 
861  j = mod(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
862 
863  ! copy the position and the error vector in the diis buffers
864 
865  IF (do_ks) THEN
866  DO ispin = 1, nspin
867  CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
868  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
869  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, qs_ot_env(ispin)%rot_mat_x)
870  END IF
871  END DO
872  END IF
873  IF (do_ener) THEN
874  DO ispin = 1, nspin
875  qs_ot_env(ispin)%ener_h_x(j, :) = qs_ot_env(ispin)%ener_x(:)
876  END DO
877  END IF
878  IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
879  qs_ot_env(1)%gnorm = 0.0_dp
880  IF (do_ks) THEN
881  DO ispin = 1, nspin
882  CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
883  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix)
884  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
885  tmp)
886  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
887  END DO
888  IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
889  WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
890  END IF
891  DO ispin = 1, nspin
892  CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -qs_ot_env(1)%ds_min)
893  END DO
894  IF (qs_ot_env(1)%settings%do_rotation) THEN
895  DO ispin = 1, nspin
896  CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, qs_ot_env(ispin)%rot_mat_gx)
897  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, &
898  tmp)
899  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
900  END DO
901  DO ispin = 1, nspin
902  CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, -qs_ot_env(1)%ds_min)
903  END DO
904  END IF
905  END IF
906  IF (do_ener) THEN
907  DO ispin = 1, nspin
908  qs_ot_env(ispin)%ener_h_e(j, :) = qs_ot_env(ispin)%ener_gx(:)
909  tmp = dot_product(qs_ot_env(ispin)%ener_h_e(j, :), qs_ot_env(ispin)%ener_gx(:))
910  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
911  qs_ot_env(ispin)%ener_h_e(j, :) = -qs_ot_env(1)%ds_min*qs_ot_env(ispin)%ener_h_e(j, :)
912  END DO
913  END IF
914  ELSE
915  qs_ot_env(1)%gnorm = 0.0_dp
916  IF (do_ks) THEN
917  DO ispin = 1, nspin
918  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
919  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
920  CALL dbcsr_add(qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
921  qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-qs_ot_env(1)%ds_min)
922  END DO
923  IF (qs_ot_env(1)%settings%do_rotation) THEN
924  DO ispin = 1, nspin
925  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
926  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
927  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, &
928  qs_ot_env(ispin)%rot_mat_gx, alpha_scalar=0.0_dp, beta_scalar=-qs_ot_env(1)%ds_min)
929  END DO
930  END IF
931  END IF
932  IF (do_ener) THEN
933  DO ispin = 1, nspin
934  tmp = dot_product(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_gx(:))
935  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
936  qs_ot_env(ispin)%ener_h_e(j, :) = -qs_ot_env(1)%ds_min*qs_ot_env(ispin)%ener_gx(:)
937  END DO
938  END IF
939  END IF
940  k = 0
941  n = 0
942  nener = 0
943  IF (do_ks) THEN
944  CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
945  DO ispin = 1, nspin
946  CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
947  k = k + itmp
948  END DO
949  END IF
950  IF (do_ener) THEN
951  DO ispin = 1, nspin
952  nener = nener + SIZE(qs_ot_env(ispin)%ener_x)
953  END DO
954  END IF
955  ! Handling the case of no free variables to optimize
956  IF (int(n, kind=int_8)*int(k, kind=int_8) + nener /= 0) THEN
957  qs_ot_env(1)%delta = sqrt(abs(qs_ot_env(1)%gnorm)/(int(n, kind=int_8)*int(k, kind=int_8) + nener))
958  qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
959  ELSE
960  qs_ot_env(1)%delta = 0.0_dp
961  qs_ot_env(1)%gradient = 0.0_dp
962  END IF
963 
964  ! make the diis matrix and solve it
965  DO i = 1, diis_bound
966  ! I think there are two possible options, with and without preconditioner
967  ! as a metric
968  ! the second option seems most logical to me, and it seems marginally faster
969  ! in some of the tests
970  IF (.false.) THEN
971  qs_ot_env(1)%ls_diis(i, j) = 0.0_dp
972  IF (do_ks) THEN
973  DO ispin = 1, nspin
974  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
975  qs_ot_env(ispin)%matrix_h_e(i)%matrix, &
976  tmp)
977  qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) + tmp
978  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
979  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, &
980  qs_ot_env(ispin)%rot_mat_h_e(i)%matrix, &
981  tmp)
982  qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) + 0.5_dp*tmp
983  END IF
984  END DO
985  END IF
986  IF (do_ener) THEN
987  DO ispin = 1, nspin
988  tmp = dot_product(qs_ot_env(ispin)%ener_h_e(j, :), qs_ot_env(ispin)%ener_h_e(i, :))
989  qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) + tmp
990  END DO
991  END IF
992  ELSE
993  qs_ot_env(1)%ls_diis(i, j) = 0.0_dp
994  IF (do_ks) THEN
995  DO ispin = 1, nspin
996  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
997  qs_ot_env(ispin)%matrix_h_e(i)%matrix, &
998  tmp)
999  qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*tmp
1000  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1001  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, &
1002  qs_ot_env(ispin)%rot_mat_h_e(i)%matrix, &
1003  tmp)
1004  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
1005  END IF
1006  END DO
1007  END IF
1008  IF (do_ener) THEN
1009  DO ispin = 1, nspin
1010  tmp = dot_product(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_h_e(i, :))
1011  qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*tmp
1012  END DO
1013  END IF
1014  END IF
1015  qs_ot_env(1)%ls_diis(j, i) = qs_ot_env(1)%ls_diis(i, j)
1016  qs_ot_env(1)%ls_diis(i, diis_bound + 1) = 1.0_dp
1017  qs_ot_env(1)%ls_diis(diis_bound + 1, i) = 1.0_dp
1018  qs_ot_env(1)%c_diis(i) = 0.0_dp
1019  END DO
1020  qs_ot_env(1)%ls_diis(diis_bound + 1, diis_bound + 1) = 0.0_dp
1021  qs_ot_env(1)%c_diis(diis_bound + 1) = 1.0_dp
1022  ! put in buffer, dgesv destroys
1023  qs_ot_env(1)%lss_diis = qs_ot_env(1)%ls_diis
1024 
1025  CALL dgesv(diis_bound + 1, 1, qs_ot_env(1)%lss_diis, diis_m + 1, qs_ot_env(1)%ipivot, &
1026  qs_ot_env(1)%c_diis, diis_m + 1, info)
1027 
1028  IF (info .NE. 0) THEN
1029  do_ot_sd = .true.
1030  WRITE (cp_logger_get_default_unit_nr(logger), *) "Singular DIIS matrix"
1031  ELSE
1032  do_ot_sd = .false.
1033  IF (do_ks) THEN
1034  DO ispin = 1, nspin
1035  ! OK, add the vectors now
1036  CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1037  DO i = 1, diis_bound
1038  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1039  qs_ot_env(ispin)%matrix_h_e(i)%matrix, &
1040  alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1041  END DO
1042  DO i = 1, diis_bound
1043  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1044  qs_ot_env(ispin)%matrix_h_x(i)%matrix, &
1045  alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1046  END DO
1047  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1048  CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
1049  DO i = 1, diis_bound
1050  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, &
1051  qs_ot_env(ispin)%rot_mat_h_e(i)%matrix, &
1052  alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1053  END DO
1054  DO i = 1, diis_bound
1055  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, &
1056  qs_ot_env(ispin)%rot_mat_h_x(i)%matrix, &
1057  alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1058  END DO
1059  END IF
1060  END DO
1061  END IF
1062  IF (do_ener) THEN
1063  DO ispin = 1, nspin
1064  qs_ot_env(ispin)%ener_x(:) = 0.0_dp
1065  DO i = 1, diis_bound
1066  qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1067  + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_e(i, :)
1068  END DO
1069  DO i = 1, diis_bound
1070  qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1071  + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_x(i, :)
1072  END DO
1073  END DO
1074  END IF
1075  qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1076  IF (qs_ot_env(1)%settings%safer_diis) THEN
1077  ! now, final check, is the step in fact in the direction of the -gradient ?
1078  ! if not we're walking towards a sadle point, and should avoid that
1079  ! the direction of the step is x_new-x_old
1080  tr_xold_gx = 0.0_dp
1081  tr_xnew_gx = 0.0_dp
1082  IF (do_ks) THEN
1083  DO ispin = 1, nspin
1084  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_x(j)%matrix, &
1085  qs_ot_env(ispin)%matrix_gx, tmp)
1086  tr_xold_gx = tr_xold_gx + tmp
1087  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, &
1088  qs_ot_env(ispin)%matrix_gx, tmp)
1089  tr_xnew_gx = tr_xnew_gx + tmp
1090  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1091  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, &
1092  qs_ot_env(ispin)%rot_mat_gx, tmp)
1093  tr_xold_gx = tr_xold_gx + 0.5_dp*tmp
1094  CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_x, &
1095  qs_ot_env(ispin)%rot_mat_gx, tmp)
1096  tr_xnew_gx = tr_xnew_gx + 0.5_dp*tmp
1097  END IF
1098  END DO
1099  END IF
1100  IF (do_ener) THEN
1101  DO ispin = 1, nspin
1102  tmp = dot_product(qs_ot_env(ispin)%ener_h_x(j, :), qs_ot_env(ispin)%ener_gx(:))
1103  tr_xold_gx = tr_xold_gx + tmp
1104  tmp = dot_product(qs_ot_env(ispin)%ener_x(:), qs_ot_env(ispin)%ener_gx(:))
1105  tr_xnew_gx = tr_xnew_gx + tmp
1106  END DO
1107  END IF
1108  overlap = (tr_xnew_gx - tr_xold_gx)
1109  ! OK, bad luck, take a SD step along the preconditioned gradient
1110  IF (overlap .GT. 0.0_dp) THEN
1111  do_ot_sd = .true.
1112  END IF
1113  END IF
1114  END IF
1115 
1116  IF (do_ot_sd) THEN
1117  qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
1118  IF (do_ks) THEN
1119  DO ispin = 1, nspin
1120  CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1121  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1122  qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
1123  1.0_dp, 1.0_dp)
1124  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1125  qs_ot_env(ispin)%matrix_h_x(j)%matrix, &
1126  1.0_dp, 1.0_dp)
1127  IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1128  CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
1129  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, &
1130  qs_ot_env(ispin)%rot_mat_h_e(j)%matrix, &
1131  1.0_dp, 1.0_dp)
1132  CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_x, &
1133  qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, &
1134  1.0_dp, 1.0_dp)
1135  END IF
1136  END DO
1137  END IF
1138  IF (do_ener) THEN
1139  DO ispin = 1, nspin
1140  qs_ot_env(ispin)%ener_x(:) = 0._dp
1141  qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_e(j, :)
1142  qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_x(j, :)
1143  END DO
1144  END IF
1145  END IF
1146 
1147  CALL timestop(handle)
1148 
1149  END SUBROUTINE ot_diis_step
1150 
1151 ! **************************************************************************************************
1152 !> \brief Energy minimizer by Broyden's method
1153 !> \param qs_ot_env variable to control minimizer behaviour
1154 !> \author Kurt Baarman (09.2010)
1155 ! **************************************************************************************************
1156  SUBROUTINE ot_broyden_step(qs_ot_env)
1157  TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
1158 
1159  INTEGER :: diis_bound, diis_m, i, ispin, itmp, j, &
1160  k, n, nspin
1161  INTEGER, ALLOCATABLE, DIMENSION(:) :: circ_index
1162  LOGICAL :: adaptive_sigma, do_ener, do_ks, &
1163  enable_flip, forget_history
1164  REAL(kind=dp) :: beta, eta, gamma, omega, sigma, &
1165  sigma_dec, sigma_min, tmp, tmp2
1166  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f, x
1167  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: g, s
1168  TYPE(cp_logger_type), POINTER :: logger
1169 
1170  eta = qs_ot_env(1)%settings%broyden_eta
1171  omega = qs_ot_env(1)%settings%broyden_omega
1172  sigma_dec = qs_ot_env(1)%settings%broyden_sigma_decrease
1173  sigma_min = qs_ot_env(1)%settings%broyden_sigma_min
1174  forget_history = qs_ot_env(1)%settings%broyden_forget_history
1175  adaptive_sigma = qs_ot_env(1)%settings%broyden_adaptive_sigma
1176  enable_flip = qs_ot_env(1)%settings%broyden_enable_flip
1177 
1178  do_ks = qs_ot_env(1)%settings%ks
1179  do_ener = qs_ot_env(1)%settings%do_ener
1180 
1181  beta = qs_ot_env(1)%settings%broyden_beta
1182  gamma = qs_ot_env(1)%settings%broyden_gamma
1183  IF (adaptive_sigma) THEN
1184  IF (qs_ot_env(1)%broyden_adaptive_sigma .LT. 0.0_dp) THEN
1185  sigma = qs_ot_env(1)%settings%broyden_sigma
1186  ELSE
1187  sigma = qs_ot_env(1)%broyden_adaptive_sigma
1188  END IF
1189  ELSE
1190  sigma = qs_ot_env(1)%settings%broyden_sigma
1191  END IF
1192 
1193  ! simplify our life....
1194  IF (do_ener .OR. .NOT. do_ks .OR. qs_ot_env(1)%settings%do_rotation) THEN
1195  cpabort("Not yet implemented")
1196  END IF
1197  !
1198  nspin = SIZE(qs_ot_env)
1199 
1200  diis_m = qs_ot_env(1)%settings%diis_m
1201 
1202  IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
1203  diis_bound = qs_ot_env(1)%diis_iter + 1
1204  ELSE
1205  diis_bound = diis_m
1206  END IF
1207 
1208  ! We want x:s, f:s and one random vector
1209  k = 2*diis_bound + 1
1210  ALLOCATE (s(k, k))
1211  ALLOCATE (g(k, k))
1212  ALLOCATE (f(k))
1213  ALLOCATE (x(k))
1214  ALLOCATE (circ_index(diis_bound))
1215  g = 0.0
1216  DO i = 1, k
1217  g(i, i) = sigma
1218  END DO
1219  s = 0.0
1220 
1221  j = mod(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
1222 
1223  DO ispin = 1, nspin
1224  CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
1225  END DO
1226 
1227  IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
1228  qs_ot_env(1)%gnorm = 0.0_dp
1229  DO ispin = 1, nspin
1230  CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
1231  qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix)
1232  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
1233  tmp)
1234  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1235  END DO
1236  IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
1237  WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1238  END IF
1239  DO ispin = 1, nspin
1240  CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -1.0_dp)
1241  END DO
1242  ELSE
1243  qs_ot_env(1)%gnorm = 0.0_dp
1244  DO ispin = 1, nspin
1245  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
1246  qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1247  CALL dbcsr_add(qs_ot_env(ispin)%matrix_h_e(j)%matrix, &
1248  qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-1.0_dp)
1249  END DO
1250  END IF
1251 
1252  k = 0
1253  n = 0
1254  CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
1255  DO ispin = 1, nspin
1256  CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
1257  k = k + itmp
1258  END DO
1259 
1260  ! Handling the case of no free variables to optimize
1261  IF (int(n, kind=int_8)*int(k, kind=int_8) /= 0) THEN
1262  qs_ot_env(1)%delta = sqrt(abs(qs_ot_env(1)%gnorm)/(int(n, kind=int_8)*int(k, kind=int_8)))
1263  qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
1264  ELSE
1265  qs_ot_env(1)%delta = 0.0_dp
1266  qs_ot_env(1)%gradient = 0.0_dp
1267  END IF
1268 
1269  IF (diis_bound == diis_m) THEN
1270  DO i = 1, diis_bound
1271  circ_index(i) = mod(j + i - 1, diis_m) + 1
1272  END DO
1273  ELSE
1274  DO i = 1, diis_bound
1275  circ_index(i) = i
1276  END DO
1277  END IF
1278 
1279  s = 0.0_dp
1280  DO ispin = 1, nspin
1281  CALL dbcsr_init_random(qs_ot_env(ispin)%matrix_x)
1282  DO i = 1, diis_bound
1283  CALL dbcsr_dot( &
1284  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1285  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, tmp)
1286  s(i, i) = s(i, i) + tmp
1287  CALL dbcsr_dot( &
1288  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1289  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, tmp)
1290  s(i + diis_bound, i + diis_bound) = s(i + diis_bound, i + diis_bound) + tmp
1291  CALL dbcsr_dot( &
1292  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1293  qs_ot_env(ispin)%matrix_x, tmp)
1294  s(i, 2*diis_bound + 1) = s(i, 2*diis_bound + 1) + tmp
1295  s(i, 2*diis_bound + 1) = s(2*diis_bound + 1, i)
1296  CALL dbcsr_dot( &
1297  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1298  qs_ot_env(ispin)%matrix_x, tmp)
1299  s(i + diis_bound, 2*diis_bound + 1) = s(i + diis_bound, 2*diis_bound + 1) + tmp
1300  s(i + diis_bound, 2*diis_bound + 1) = s(2*diis_bound + 1, diis_bound + i)
1301  DO k = (i + 1), diis_bound
1302  CALL dbcsr_dot( &
1303  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1304  qs_ot_env(ispin)%matrix_h_x(circ_index(k))%matrix, &
1305  tmp)
1306  s(i, k) = s(i, k) + tmp
1307  s(k, i) = s(i, k)
1308  CALL dbcsr_dot( &
1309  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1310  qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix, &
1311  tmp)
1312  s(diis_bound + i, diis_bound + k) = s(diis_bound + i, diis_bound + k) + tmp
1313  s(diis_bound + k, diis_bound + i) = s(diis_bound + i, diis_bound + k)
1314  END DO
1315  DO k = 1, diis_bound
1316  CALL dbcsr_dot( &
1317  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1318  qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix, tmp)
1319  s(i, k + diis_bound) = s(i, k + diis_bound) + tmp
1320  s(k + diis_bound, i) = s(i, k + diis_bound)
1321  END DO
1322  END DO
1323  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_x, tmp)
1324  s(2*diis_bound + 1, 2*diis_bound + 1) = s(2*diis_bound + 1, 2*diis_bound + 1) + tmp
1325  END DO
1326 
1327  ! normalize
1328  k = 2*diis_bound + 1
1329  tmp = sqrt(s(k, k))
1330  s(k, :) = s(k, :)/tmp
1331  s(:, k) = s(:, k)/tmp
1332 
1333  IF (diis_bound .GT. 1) THEN
1334  tmp = 0.0_dp
1335  tmp2 = 0.0_dp
1336  i = diis_bound
1337  DO ispin = 1, nspin
1338  ! dot product of differences
1339  CALL dbcsr_dot( &
1340  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1341  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1342  tmp)
1343  tmp2 = tmp2 + tmp
1344  CALL dbcsr_dot( &
1345  qs_ot_env(ispin)%matrix_h_x(circ_index(i - 1))%matrix, &
1346  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1347  tmp)
1348  tmp2 = tmp2 - tmp
1349  CALL dbcsr_dot( &
1350  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1351  qs_ot_env(ispin)%matrix_h_e(circ_index(i - 1))%matrix, &
1352  tmp)
1353  tmp2 = tmp2 - tmp
1354  CALL dbcsr_dot( &
1355  qs_ot_env(ispin)%matrix_h_x(circ_index(i - 1))%matrix, &
1356  qs_ot_env(ispin)%matrix_h_e(circ_index(i - 1))%matrix, &
1357  tmp)
1358  tmp2 = tmp2 + tmp
1359  END DO
1360  qs_ot_env(1)%c_broy(i - 1) = tmp2
1361  END IF
1362 
1363  qs_ot_env(1)%energy_h(j) = qs_ot_env(1)%etotal
1364 
1365  ! If we went uphill, do backtracking line search
1366  i = minloc(qs_ot_env(1)%energy_h(1:diis_bound), dim=1)
1367  IF (i .NE. j) THEN
1368  sigma = sigma_dec*sigma
1369  qs_ot_env(1)%OT_METHOD_FULL = "OT BTRK"
1370  DO ispin = 1, nspin
1371  CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1372  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1373  qs_ot_env(ispin)%matrix_h_x(i)%matrix, &
1374  alpha_scalar=1.0_dp, beta_scalar=(1.0_dp - gamma))
1375  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1376  qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix, &
1377  alpha_scalar=1.0_dp, beta_scalar=gamma)
1378  END DO
1379  ELSE
1380  ! Construct G
1381  DO i = 2, diis_bound
1382  f = 0.0
1383  x = 0.0
1384  ! f is df_i
1385  x(i) = 1.0
1386  x(i - 1) = -1.0
1387  ! x is dx_i
1388  f(diis_bound + i) = 1.0
1389  f(diis_bound + i - 1) = -1.0
1390  tmp = 1.0_dp
1391  ! We want a pos def Hessian
1392  IF (enable_flip) THEN
1393  IF (qs_ot_env(1)%c_broy(i - 1) .GT. 0) THEN
1394  !qs_ot_env(1)%OT_METHOD_FULL="OT FLIP"
1395  tmp = -1.0_dp
1396  END IF
1397  END IF
1398 
1399  ! get dx-Gdf
1400  x(:) = tmp*x - matmul(g, f)
1401  ! dfSdf
1402  ! we calculate matmul(S, f) twice. They're small...
1403  tmp = dot_product(f, matmul(s, f))
1404  ! NOTE THAT S IS SYMMETRIC !!!
1405  f(:) = matmul(s, f)/tmp
1406  ! the spread is an outer vector product
1407  g(:, :) = g + spread(x, dim=2, ncopies=SIZE(f))*spread(f, dim=1, ncopies=SIZE(x))
1408  END DO
1409  f = 0.0_dp
1410  f(2*diis_bound) = 1.0_dp
1411  x(:) = -beta*matmul(g, f)
1412 
1413  ! OK, add the vectors now, this sums up to the proposed step
1414  DO ispin = 1, nspin
1415  CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1416  DO i = 1, diis_bound
1417  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1418  qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1419  alpha_scalar=1.0_dp, beta_scalar=-x(i + diis_bound))
1420  END DO
1421  DO i = 1, diis_bound
1422  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1423  qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1424  alpha_scalar=1.0_dp, beta_scalar=x(i))
1425  END DO
1426  END DO
1427 
1428  IF (adaptive_sigma) THEN
1429  tmp = new_sigma(g, s, diis_bound)
1430  !tmp = tmp * qs_ot_env ( 1 ) % settings % broyden_sigma
1431  tmp = tmp*eta
1432  sigma = min(omega*sigma, tmp)
1433  END IF
1434 
1435  ! compute the inner product of direction of the step and gradient
1436  tmp = 0.0_dp
1437  DO ispin = 1, nspin
1438  CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
1439  qs_ot_env(ispin)%matrix_x, &
1440  tmp2)
1441  tmp = tmp + tmp2
1442  END DO
1443 
1444  DO ispin = 1, nspin
1445  ! if the direction of the step is not in direction of the gradient,
1446  ! change step sign
1447  IF (tmp .GE. 0.0_dp) THEN
1448  qs_ot_env(1)%OT_METHOD_FULL = "OT TURN"
1449  IF (forget_history) THEN
1450  qs_ot_env(1)%diis_iter = 0
1451  END IF
1452  sigma = sigma*sigma_dec
1453  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1454  qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix, &
1455  alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
1456  ELSE
1457  CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, &
1458  qs_ot_env(ispin)%matrix_h_x(circ_index(diis_bound))%matrix, &
1459  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1460  END IF
1461  END DO
1462  END IF
1463 
1464  ! get rid of S, G, f, x, circ_index for next round
1465  DEALLOCATE (s, g, f, x, circ_index)
1466 
1467  ! update for next round
1468  qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1469  qs_ot_env(1)%broyden_adaptive_sigma = max(sigma, sigma_min)
1470 
1471  END SUBROUTINE ot_broyden_step
1472 
1473 ! **************************************************************************************************
1474 !> \brief ...
1475 !> \param G ...
1476 !> \param S ...
1477 !> \param n ...
1478 !> \return ...
1479 ! **************************************************************************************************
1480  FUNCTION new_sigma(G, S, n) RESULT(sigma)
1481 !
1482 ! Calculate new sigma from eigenvalues of full size G by Arnoldi.
1483 !
1484 ! **************************************************************************************************
1485 
1486  REAL(kind=dp), DIMENSION(:, :) :: g, s
1487  INTEGER :: n
1488  REAL(kind=dp) :: sigma
1489 
1490  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigv
1491  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: h
1492 
1493  ALLOCATE (h(n, n))
1494  CALL hess_g(g, s, h, n)
1495  ALLOCATE (eigv(n))
1496  CALL diamat_all(h(1:n, 1:n), eigv)
1497 
1498  SELECT CASE (1)
1499  CASE (1)
1500  ! This estimator seems to work well. No theory.
1501  sigma = sum(abs(eigv**2))/sum(abs(eigv))
1502  CASE (2)
1503  ! Estimator based on Frobenius norm minimizer
1504  sigma = sum(abs(eigv))/max(1, SIZE(eigv))
1505  CASE (3)
1506  ! Estimator based on induced 2-norm
1507  sigma = (maxval(abs(eigv)) + minval(abs(eigv)))*0.5_dp
1508  END SELECT
1509 
1510  DEALLOCATE (h, eigv)
1511  END FUNCTION new_sigma
1512 
1513 ! **************************************************************************************************
1514 !> \brief ...
1515 !> \param G ...
1516 !> \param S ...
1517 !> \param H ...
1518 !> \param n ...
1519 ! **************************************************************************************************
1520  SUBROUTINE hess_g(G, S, H, n)
1521 !
1522 ! Make a hessenberg out of G into H. Cf Arnoldi.
1523 ! Inner product is weighted by S.
1524 ! Possible lucky breakdown at n.
1525 !
1526 ! **************************************************************************************************
1527  REAL(kind=dp), DIMENSION(:, :) :: g, s, h
1528  INTEGER :: n
1529 
1530  INTEGER :: i, j, k
1531  REAL(kind=dp) :: tmp
1532  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: v
1533  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: q
1534 
1535  i = SIZE(g, 1)
1536  k = SIZE(h, 1)
1537  ALLOCATE (q(i, k))
1538  ALLOCATE (v(i))
1539  h = 0.0_dp
1540  q = 0.0_dp
1541 
1542  q(:, 1) = 1.0_dp
1543  tmp = sqrt(dot_product(q(:, 1), matmul(s, q(:, 1))))
1544  q(:, :) = q(:, :)/tmp
1545 
1546  DO i = 1, k
1547  v(:) = matmul(g, q(:, i))
1548  DO j = 1, i
1549  h(j, i) = dot_product(q(:, j), matmul(s, v))
1550  v(:) = v - h(j, i)*q(:, j)
1551  END DO
1552  IF (i .LT. k) THEN
1553  tmp = dot_product(v, matmul(s, v))
1554  IF (tmp .LE. 0.0_dp) THEN
1555  n = i
1556  EXIT
1557  END IF
1558  tmp = sqrt(tmp)
1559  ! Lucky breakdown
1560  IF (abs(tmp) .LT. 1e-9_dp) THEN
1561  n = i
1562  EXIT
1563  END IF
1564  h(i + 1, i) = tmp
1565  q(:, i + 1) = v/h(i + 1, i)
1566  END IF
1567  END DO
1568 
1569  DEALLOCATE (q, v)
1570  END SUBROUTINE hess_g
1571 
1572 END MODULE qs_ot_minimizer
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:567
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
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:372
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::ot_mini
subroutine, public ot_mini(qs_ot_env, matrix_hc)
...
Definition: qs_ot_minimizer.F:63
qs_ot_types
orbital transformations
Definition: qs_ot_types.F:15
qs_ot
orbital transformations
Definition: qs_ot.F:15
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:704
qs_ot::qs_ot_get_derivative
subroutine, public qs_ot_get_derivative(matrix_hc, matrix_x, matrix_sx, matrix_gx, qs_ot_env)
...
Definition: qs_ot.F:1074
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17