Line data Source code
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 78554 : 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 78554 : CALL timeset(routineN, handle)
73 :
74 78554 : nspin = SIZE(qs_ot_env)
75 :
76 78554 : do_ks = qs_ot_env(1)%settings%ks
77 78554 : do_ener = qs_ot_env(1)%settings%do_ener
78 :
79 78554 : qs_ot_env(1)%OT_METHOD_FULL = ""
80 :
81 : ! compute the gradient for the variables x
82 78554 : IF (.NOT. qs_ot_env(1)%energy_only) THEN
83 61135 : qs_ot_env(1)%gradient = 0.0_dp
84 129099 : DO ispin = 1, nspin
85 67964 : IF (do_ks) THEN
86 133076 : 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 65112 : 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 2852 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
95 : CASE DEFAULT
96 67964 : CPABORT("ALGORITHM NYI")
97 : END SELECT
98 : END IF
99 : ! and also the gradient along the direction
100 129099 : IF (qs_ot_env(1)%use_dx) THEN
101 61440 : IF (do_ks) THEN
102 61440 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
103 61440 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
104 61440 : IF (qs_ot_env(1)%settings%do_rotation) THEN
105 1404 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
106 1404 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + 0.5_dp*tmp
107 : END IF
108 : END IF
109 61440 : IF (do_ener) THEN
110 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
111 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
112 : END IF
113 : ELSE
114 6524 : IF (do_ks) THEN
115 6524 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
116 6524 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
117 6524 : IF (qs_ot_env(1)%settings%do_rotation) THEN
118 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
119 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - 0.5_dp*tmp
120 : END IF
121 : END IF
122 6524 : IF (do_ener) THEN
123 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
124 0 : 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 114770 : SELECT CASE (qs_ot_env(1)%settings%OT_METHOD)
131 : CASE ("CG")
132 36216 : IF (current_point_is_fine(qs_ot_env)) THEN
133 18815 : qs_ot_env(1)%OT_METHOD_FULL = "OT CG"
134 18815 : CALL ot_new_cg_direction(qs_ot_env)
135 18815 : qs_ot_env(1)%line_search_count = 0
136 : ELSE
137 17401 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
138 : END IF
139 36216 : CALL do_line_search(qs_ot_env)
140 : CASE ("SD")
141 42 : IF (current_point_is_fine(qs_ot_env)) THEN
142 24 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
143 24 : CALL ot_new_sd_direction(qs_ot_env)
144 24 : qs_ot_env(1)%line_search_count = 0
145 : ELSE
146 18 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
147 : END IF
148 42 : CALL do_line_search(qs_ot_env)
149 : CASE ("DIIS")
150 42118 : qs_ot_env(1)%OT_METHOD_FULL = "OT DIIS"
151 42118 : CALL ot_diis_step(qs_ot_env)
152 : CASE ("BROY")
153 178 : qs_ot_env(1)%OT_METHOD_FULL = "OT BROY"
154 178 : CALL ot_broyden_step(qs_ot_env)
155 : CASE DEFAULT
156 78554 : CPABORT("OT_METHOD NYI")
157 : END SELECT
158 :
159 78554 : CALL timestop(handle)
160 :
161 78554 : 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 36258 : 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 36258 : res = .FALSE.
177 :
178 : ! only if we have a gradient it can be fine
179 36258 : IF (.NOT. qs_ot_env(1)%energy_only) THEN
180 :
181 : ! we have not yet started with the line search
182 18839 : IF (qs_ot_env(1)%line_search_count .EQ. 0) THEN
183 36258 : res = .TRUE.
184 : RETURN
185 : END IF
186 :
187 15888 : IF (qs_ot_env(1)%line_search_might_be_done) THEN
188 : ! here we put the more complicated logic later
189 15888 : res = .TRUE.
190 15888 : 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 36258 : SUBROUTINE do_line_search(qs_ot_env)
205 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
206 :
207 36642 : SELECT CASE (qs_ot_env(1)%settings%line_search_method)
208 : CASE ("GOLD")
209 384 : CALL do_line_search_gold(qs_ot_env)
210 : CASE ("3PNT")
211 1186 : CALL do_line_search_3pnt(qs_ot_env)
212 : CASE ("2PNT")
213 34678 : CALL do_line_search_2pnt(qs_ot_env)
214 : CASE ("NONE")
215 10 : CALL do_line_search_none(qs_ot_env)
216 : CASE DEFAULT
217 36258 : CPABORT("NYI")
218 : END SELECT
219 36258 : 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 36258 : 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 36258 : CALL timeset(routineN, handle)
238 :
239 36258 : nspin = SIZE(qs_ot_env)
240 :
241 36258 : do_ks = qs_ot_env(1)%settings%ks
242 36258 : 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 36258 : IF (qs_ot_env(1)%use_dx) THEN
247 36246 : IF (do_ks) THEN
248 79384 : DO ispin = 1, nspin
249 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_dx, &
250 43138 : alpha_scalar=1.0_dp, beta_scalar=ds)
251 79384 : 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 2034 : alpha_scalar=1.0_dp, beta_scalar=ds)
254 : END IF
255 : END DO
256 : END IF
257 36246 : IF (do_ener) THEN
258 0 : DO ispin = 1, nspin
259 0 : 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 12 : IF (do_ks) THEN
264 24 : DO ispin = 1, nspin
265 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_gx, &
266 12 : alpha_scalar=1.0_dp, beta_scalar=-ds)
267 24 : 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 0 : alpha_scalar=1.0_dp, beta_scalar=-ds)
270 : END IF
271 : END DO
272 : END IF
273 12 : IF (do_ener) THEN
274 0 : DO ispin = 1, nspin
275 0 : 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 36258 : CALL timestop(handle)
280 36258 : 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 384 : 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 384 : CALL timeset(routineN, handle)
298 :
299 384 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
300 384 : count = qs_ot_env(1)%line_search_count
301 384 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
302 384 : qs_ot_env(1)%energy_only = .TRUE.
303 :
304 384 : 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 0 : CPABORT("MAX ITER EXCEEDED : FATAL")
308 : END IF
309 :
310 384 : IF (qs_ot_env(1)%line_search_count .EQ. 1) THEN
311 30 : qs_ot_env(1)%line_search_left = 1
312 30 : qs_ot_env(1)%line_search_right = 0
313 30 : qs_ot_env(1)%line_search_mid = 1
314 30 : qs_ot_env(1)%ot_pos(1) = 0.0_dp
315 30 : qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
316 30 : qs_ot_env(1)%ot_pos(2) = qs_ot_env(1)%ds_min/gold_sec
317 : ELSE
318 354 : 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 354 : IF (qs_ot_env(1)%line_search_right .EQ. 0) THEN ! we do not yet have the right bracket
323 40 : IF (qs_ot_env(1)%ot_energy(count - 1) .LT. qs_ot_env(1)%ot_energy(count)) THEN
324 26 : 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 26 : qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))*gold_sec
328 : ELSE
329 14 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
330 14 : qs_ot_env(1)%line_search_mid = count
331 14 : 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 314 : IF (qs_ot_env(1)%ot_pos(count) .LT. qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
336 140 : IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
337 42 : qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
338 42 : qs_ot_env(1)%line_search_mid = count
339 : ELSE
340 98 : qs_ot_env(1)%line_search_left = count
341 : END IF
342 : ELSE
343 174 : IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
344 66 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
345 66 : qs_ot_env(1)%line_search_mid = count
346 : ELSE
347 108 : 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 314 : 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 154 : - 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 160 : - 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 314 : 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 20 : qs_ot_env(1)%energy_only = .FALSE.
373 20 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
374 : END IF
375 : END IF
376 : END IF
377 384 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
378 384 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
379 :
380 384 : CALL take_step(ds, qs_ot_env)
381 :
382 384 : CALL timestop(handle)
383 :
384 384 : END SUBROUTINE do_line_search_gold
385 :
386 : ! **************************************************************************************************
387 : !> \brief ...
388 : !> \param qs_ot_env ...
389 : ! **************************************************************************************************
390 1186 : 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 1186 : CALL timeset(routineN, handle)
401 :
402 1186 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
403 1186 : qs_ot_env(1)%energy_only = .TRUE.
404 :
405 : ! a three point interpolation based on the energy
406 1186 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
407 1186 : count = qs_ot_env(1)%line_search_count
408 1186 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
409 410 : SELECT CASE (count)
410 : CASE (1)
411 410 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
412 410 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*0.8_dp
413 : CASE (2)
414 394 : IF (qs_ot_env(1)%OT_energy(count) .GT. qs_ot_env(1)%OT_energy(count - 1)) THEN
415 38 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*0.5_dp
416 : ELSE
417 356 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*1.4_dp
418 : END IF
419 : CASE (3)
420 382 : xa = qs_ot_env(1)%OT_pos(1)
421 382 : xb = qs_ot_env(1)%OT_pos(2)
422 382 : xc = qs_ot_env(1)%OT_pos(3)
423 382 : fa = qs_ot_env(1)%OT_energy(1)
424 382 : fb = qs_ot_env(1)%OT_energy(2)
425 382 : fc = qs_ot_env(1)%OT_energy(3)
426 382 : nom = (xb - xa)**2*(fb - fc) - (xb - xc)**2*(fb - fa)
427 382 : denom = (xb - xa)*(fb - fc) - (xb - xc)*(fb - fa)
428 382 : IF (ABS(denom) .LE. 1.0E-18_dp*MAX(ABS(fb - fc), ABS(fb - fa))) THEN
429 : pos = xb
430 : ELSE
431 382 : 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 382 : (pos - xc)*(pos - xa)*fb/((xb - xc)*(xb - xa))
436 382 : 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 1900 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:3))*4.0_dp))
440 : ELSE ! just take an extended step
441 10 : qs_ot_env(1)%OT_pos(count + 1) = MAXVAL(qs_ot_env(1)%OT_pos(1:3))*2.0_dp
442 : END IF
443 382 : qs_ot_env(1)%energy_only = .FALSE.
444 382 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
445 : CASE DEFAULT
446 1186 : CPABORT("NYI")
447 : END SELECT
448 1186 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
449 1186 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
450 :
451 1186 : CALL take_step(ds, qs_ot_env)
452 :
453 1186 : CALL timestop(handle)
454 :
455 1186 : END SUBROUTINE do_line_search_3pnt
456 :
457 : ! **************************************************************************************************
458 : !> \brief ...
459 : !> \param qs_ot_env ...
460 : ! **************************************************************************************************
461 34678 : 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 34678 : CALL timeset(routineN, handle)
471 :
472 34678 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
473 34678 : qs_ot_env(1)%energy_only = .TRUE.
474 :
475 : ! a three point interpolation based on the energy
476 34678 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
477 34678 : count = qs_ot_env(1)%line_search_count
478 34678 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
479 18389 : SELECT CASE (count)
480 : CASE (1)
481 18389 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
482 18389 : qs_ot_env(1)%ot_grad(count) = qs_ot_env(1)%gradient
483 18389 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*1.0_dp
484 : CASE (2)
485 16289 : x0 = 0.0_dp
486 16289 : c = qs_ot_env(1)%ot_energy(1)
487 16289 : b = qs_ot_env(1)%ot_grad(1)
488 16289 : x1 = qs_ot_env(1)%ot_pos(2)
489 16289 : a = (qs_ot_env(1)%ot_energy(2) - b*x1 - c)/(x1**2)
490 16289 : IF (a .LE. 0.0_dp) a = 1.0E-15_dp
491 16289 : pos = -b/(2.0_dp*a)
492 16289 : val = a*pos**2 + b*pos + c
493 16289 : qs_ot_env(1)%energy_only = .FALSE.
494 16289 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
495 16289 : 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 64820 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:2))*4.0_dp))
500 : ELSE ! just take an extended step
501 336 : 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 34678 : CPABORT("NYI")
505 : END SELECT
506 34678 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
507 34678 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
508 :
509 34678 : CALL take_step(ds, qs_ot_env)
510 :
511 34678 : CALL timestop(handle)
512 :
513 34678 : END SUBROUTINE do_line_search_2pnt
514 :
515 : ! **************************************************************************************************
516 : !> \brief ...
517 : !> \param qs_ot_env ...
518 : ! **************************************************************************************************
519 10 : SUBROUTINE do_line_search_none(qs_ot_env)
520 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
521 :
522 10 : CALL take_step(qs_ot_env(1)%ds_min, qs_ot_env)
523 :
524 10 : 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 24 : 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 24 : CALL timeset(routineN, handle)
546 :
547 : !***SCP
548 :
549 24 : nspin = SIZE(qs_ot_env)
550 24 : logger => cp_get_default_logger()
551 24 : do_ks = qs_ot_env(1)%settings%ks
552 24 : do_ener = qs_ot_env(1)%settings%do_ener
553 :
554 24 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
555 20 : IF (.NOT. qs_ot_env(1)%use_dx) CPABORT("use dx")
556 20 : qs_ot_env(1)%gnorm = 0.0_dp
557 20 : IF (do_ks) THEN
558 40 : DO ispin = 1, nspin
559 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
560 20 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx)
561 20 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
562 40 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
563 : END DO
564 20 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
565 0 : logger => cp_get_default_logger()
566 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
567 : END IF
568 40 : DO ispin = 1, nspin
569 40 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
570 : END DO
571 20 : IF (qs_ot_env(1)%settings%do_rotation) THEN
572 0 : DO ispin = 1, nspin
573 : ! right now no preconditioner yet
574 0 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
575 0 : 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 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
578 : END DO
579 0 : DO ispin = 1, nspin
580 0 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
581 : END DO
582 : END IF
583 : END IF
584 20 : IF (do_ener) THEN
585 0 : DO ispin = 1, nspin
586 0 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_gx
587 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_dx, qs_ot_env(ispin)%ener_gx)
588 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
589 0 : qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_dx
590 : END DO
591 : END IF
592 : ELSE
593 4 : qs_ot_env(1)%gnorm = 0.0_dp
594 4 : IF (do_ks) THEN
595 8 : DO ispin = 1, nspin
596 4 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
597 8 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
598 : END DO
599 4 : IF (qs_ot_env(1)%settings%do_rotation) THEN
600 0 : DO ispin = 1, nspin
601 0 : 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 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
604 : END DO
605 : END IF
606 : END IF
607 4 : IF (do_ener) THEN
608 0 : DO ispin = 1, nspin
609 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
610 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
611 : END DO
612 : END IF
613 : END IF
614 :
615 24 : k = 0
616 24 : n = 0
617 24 : nener = 0
618 24 : IF (do_ks) THEN
619 24 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
620 48 : DO ispin = 1, nspin
621 24 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
622 48 : k = k + itmp
623 : END DO
624 : END IF
625 24 : IF (do_ener) THEN
626 0 : DO ispin = 1, nspin
627 0 : 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 24 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
632 22 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
633 22 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
634 : ELSE
635 2 : qs_ot_env(1)%delta = 0.0_dp
636 2 : qs_ot_env(1)%gradient = 0.0_dp
637 : END IF
638 :
639 24 : CALL timestop(handle)
640 :
641 24 : 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 18815 : 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 18815 : CALL timeset(routineN, handle)
663 :
664 18815 : nspin = SIZE(qs_ot_env)
665 18815 : logger => cp_get_default_logger()
666 :
667 18815 : do_ks = qs_ot_env(1)%settings%ks
668 18815 : do_ener = qs_ot_env(1)%settings%do_ener
669 :
670 18815 : gnorm_cross = 0.0_dp
671 18815 : IF (do_ks) THEN
672 41469 : DO ispin = 1, nspin
673 22654 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
674 41469 : gnorm_cross = gnorm_cross + tmp
675 : END DO
676 18815 : IF (qs_ot_env(1)%settings%do_rotation) THEN
677 2024 : DO ispin = 1, nspin
678 1038 : 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 2024 : gnorm_cross = gnorm_cross + 0.5_dp*tmp
681 : END DO
682 : END IF
683 : END IF
684 18815 : IF (do_ener) THEN
685 0 : DO ispin = 1, nspin
686 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
687 0 : gnorm_cross = gnorm_cross + tmp
688 : END DO
689 : END IF
690 :
691 18815 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
692 :
693 28377 : DO ispin = 1, nspin
694 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
695 28377 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
696 : END DO
697 12509 : qs_ot_env(1)%gnorm = 0.0_dp
698 12509 : IF (do_ks) THEN
699 28377 : DO ispin = 1, nspin
700 15868 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
701 28377 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
702 : END DO
703 12509 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
704 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
705 : END IF
706 28377 : DO ispin = 1, nspin
707 28377 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
708 : END DO
709 12509 : IF (qs_ot_env(1)%settings%do_rotation) THEN
710 1976 : DO ispin = 1, nspin
711 : ! right now no preconditioner yet
712 1014 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
713 1014 : 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 1976 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
716 : END DO
717 1976 : DO ispin = 1, nspin
718 1976 : 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 12509 : IF (do_ener) THEN
723 0 : DO ispin = 1, nspin
724 0 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
725 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
726 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
727 0 : qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_gx_old
728 : END DO
729 : END IF
730 : ELSE
731 6306 : IF (do_ks) THEN
732 6306 : qs_ot_env(1)%gnorm = 0.0_dp
733 13092 : DO ispin = 1, nspin
734 6786 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
735 6786 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
736 13092 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx)
737 : END DO
738 6306 : IF (qs_ot_env(1)%settings%do_rotation) THEN
739 48 : DO ispin = 1, nspin
740 24 : 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 24 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
743 48 : 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 6306 : IF (do_ener) THEN
748 0 : DO ispin = 1, nspin
749 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
750 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
751 0 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
752 : END DO
753 : END IF
754 : END IF
755 :
756 18815 : k = 0
757 18815 : n = 0
758 18815 : nener = 0
759 18815 : IF (do_ks) THEN
760 18815 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
761 41469 : DO ispin = 1, nspin
762 22654 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
763 41469 : k = k + itmp
764 : END DO
765 : END IF
766 18815 : IF (do_ener) THEN
767 0 : DO ispin = 1, nspin
768 0 : 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 18815 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
773 18681 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
774 18681 : beta_pr = (qs_ot_env(1)%gnorm - gnorm_cross)/qs_ot_env(1)%gnorm_old
775 : ELSE
776 134 : qs_ot_env(1)%delta = 0.0_dp
777 134 : beta_pr = 0.0_dp
778 : END IF
779 18815 : beta_pr = MAX(beta_pr, 0.0_dp) ! reset to SD
780 :
781 18815 : test_down = 0.0_dp
782 18815 : IF (do_ks) THEN
783 41469 : DO ispin = 1, nspin
784 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
785 22654 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
786 22654 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
787 22654 : test_down = test_down + tmp
788 41469 : 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 1038 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
791 1038 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
792 1038 : test_down = test_down + 0.5_dp*tmp
793 : END IF
794 : END DO
795 : END IF
796 18815 : IF (do_ener) THEN
797 0 : DO ispin = 1, nspin
798 0 : qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - qs_ot_env(ispin)%ener_gx
799 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
800 0 : test_down = test_down + tmp
801 : END DO
802 : END IF
803 :
804 18815 : IF (test_down .GE. 0.0_dp) THEN ! reset to SD
805 635 : beta_pr = 0.0_dp
806 635 : IF (do_ks) THEN
807 1459 : DO ispin = 1, nspin
808 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
809 824 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
810 1459 : IF (qs_ot_env(1)%settings%do_rotation) THEN
811 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
812 10 : 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 635 : IF (do_ener) THEN
817 0 : DO ispin = 1, nspin
818 0 : 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 18815 : qs_ot_env(1)%gradient = beta_pr*qs_ot_env(1)%gradient - qs_ot_env(1)%gnorm
824 18815 : qs_ot_env(1)%gnorm_old = qs_ot_env(1)%gnorm
825 :
826 18815 : CALL timestop(handle)
827 :
828 18815 : END SUBROUTINE ot_new_cg_direction
829 :
830 : ! **************************************************************************************************
831 : !> \brief ...
832 : !> \param qs_ot_env ...
833 : ! **************************************************************************************************
834 42118 : 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 42118 : CALL timeset(routineN, handle)
846 :
847 42118 : logger => cp_get_default_logger()
848 :
849 42118 : do_ks = qs_ot_env(1)%settings%ks
850 42118 : do_ener = qs_ot_env(1)%settings%do_ener
851 42118 : nspin = SIZE(qs_ot_env)
852 :
853 42118 : diis_m = qs_ot_env(1)%settings%diis_m
854 :
855 42118 : IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
856 25828 : diis_bound = qs_ot_env(1)%diis_iter + 1
857 : ELSE
858 : diis_bound = diis_m
859 : END IF
860 :
861 42118 : 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 42118 : IF (do_ks) THEN
866 87226 : DO ispin = 1, nspin
867 45108 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
868 87226 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
869 366 : 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 42118 : IF (do_ener) THEN
874 0 : DO ispin = 1, nspin
875 0 : qs_ot_env(ispin)%ener_h_x(j, :) = qs_ot_env(ispin)%ener_x(:)
876 : END DO
877 : END IF
878 42118 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
879 35598 : qs_ot_env(1)%gnorm = 0.0_dp
880 35598 : IF (do_ks) THEN
881 74186 : DO ispin = 1, nspin
882 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
883 38588 : 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 38588 : tmp)
886 74186 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
887 : END DO
888 35598 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
889 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
890 : END IF
891 74186 : DO ispin = 1, nspin
892 74186 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -qs_ot_env(1)%ds_min)
893 : END DO
894 35598 : IF (qs_ot_env(1)%settings%do_rotation) THEN
895 732 : DO ispin = 1, nspin
896 366 : 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 366 : tmp)
899 732 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
900 : END DO
901 732 : DO ispin = 1, nspin
902 732 : 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 35598 : IF (do_ener) THEN
907 0 : DO ispin = 1, nspin
908 0 : qs_ot_env(ispin)%ener_h_e(j, :) = qs_ot_env(ispin)%ener_gx(:)
909 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_e(j, :), qs_ot_env(ispin)%ener_gx(:))
910 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
911 0 : 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 6520 : qs_ot_env(1)%gnorm = 0.0_dp
916 6520 : IF (do_ks) THEN
917 13040 : DO ispin = 1, nspin
918 6520 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
919 6520 : 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 13040 : qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-qs_ot_env(1)%ds_min)
922 : END DO
923 6520 : IF (qs_ot_env(1)%settings%do_rotation) THEN
924 0 : DO ispin = 1, nspin
925 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
926 0 : 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 0 : 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 6520 : IF (do_ener) THEN
933 0 : DO ispin = 1, nspin
934 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_gx(:))
935 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
936 0 : 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 42118 : k = 0
941 42118 : n = 0
942 42118 : nener = 0
943 42118 : IF (do_ks) THEN
944 42118 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
945 87226 : DO ispin = 1, nspin
946 45108 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
947 87226 : k = k + itmp
948 : END DO
949 : END IF
950 42118 : IF (do_ener) THEN
951 0 : DO ispin = 1, nspin
952 0 : 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 42118 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
957 42092 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
958 42092 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
959 : ELSE
960 26 : qs_ot_env(1)%delta = 0.0_dp
961 26 : qs_ot_env(1)%gradient = 0.0_dp
962 : END IF
963 :
964 : ! make the diis matrix and solve it
965 256169 : 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 214051 : qs_ot_env(1)%ls_diis(i, j) = 0.0_dp
994 214051 : IF (do_ks) THEN
995 444372 : 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 230321 : tmp)
999 230321 : qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*tmp
1000 444372 : 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 1976 : tmp)
1004 1976 : 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 214051 : IF (do_ener) THEN
1009 0 : DO ispin = 1, nspin
1010 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_h_e(i, :))
1011 0 : 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 214051 : qs_ot_env(1)%ls_diis(j, i) = qs_ot_env(1)%ls_diis(i, j)
1016 214051 : qs_ot_env(1)%ls_diis(i, diis_bound + 1) = 1.0_dp
1017 214051 : qs_ot_env(1)%ls_diis(diis_bound + 1, i) = 1.0_dp
1018 256169 : qs_ot_env(1)%c_diis(i) = 0.0_dp
1019 : END DO
1020 42118 : qs_ot_env(1)%ls_diis(diis_bound + 1, diis_bound + 1) = 0.0_dp
1021 42118 : qs_ot_env(1)%c_diis(diis_bound + 1) = 1.0_dp
1022 : ! put in buffer, dgesv destroys
1023 3198114 : 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 42118 : qs_ot_env(1)%c_diis, diis_m + 1, info)
1027 :
1028 42118 : IF (info .NE. 0) THEN
1029 0 : do_ot_sd = .TRUE.
1030 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "Singular DIIS matrix"
1031 : ELSE
1032 42118 : do_ot_sd = .FALSE.
1033 42118 : IF (do_ks) THEN
1034 87226 : DO ispin = 1, nspin
1035 : ! OK, add the vectors now
1036 45108 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1037 275429 : 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 275429 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1041 : END DO
1042 275429 : 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 275429 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1046 : END DO
1047 87226 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1048 366 : CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
1049 2342 : 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 2342 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1053 : END DO
1054 2342 : 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 2342 : 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 42118 : IF (do_ener) THEN
1063 0 : DO ispin = 1, nspin
1064 0 : qs_ot_env(ispin)%ener_x(:) = 0.0_dp
1065 0 : DO i = 1, diis_bound
1066 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1067 0 : + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_e(i, :)
1068 : END DO
1069 0 : DO i = 1, diis_bound
1070 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1071 0 : + 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 42118 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1076 42118 : 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 42118 : tr_xold_gx = 0.0_dp
1081 42118 : tr_xnew_gx = 0.0_dp
1082 42118 : IF (do_ks) THEN
1083 87226 : DO ispin = 1, nspin
1084 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_x(j)%matrix, &
1085 45108 : qs_ot_env(ispin)%matrix_gx, tmp)
1086 45108 : tr_xold_gx = tr_xold_gx + tmp
1087 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, &
1088 45108 : qs_ot_env(ispin)%matrix_gx, tmp)
1089 45108 : tr_xnew_gx = tr_xnew_gx + tmp
1090 87226 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1091 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, &
1092 366 : qs_ot_env(ispin)%rot_mat_gx, tmp)
1093 366 : tr_xold_gx = tr_xold_gx + 0.5_dp*tmp
1094 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_x, &
1095 366 : qs_ot_env(ispin)%rot_mat_gx, tmp)
1096 366 : tr_xnew_gx = tr_xnew_gx + 0.5_dp*tmp
1097 : END IF
1098 : END DO
1099 : END IF
1100 42118 : IF (do_ener) THEN
1101 0 : DO ispin = 1, nspin
1102 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_x(j, :), qs_ot_env(ispin)%ener_gx(:))
1103 0 : tr_xold_gx = tr_xold_gx + tmp
1104 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_x(:), qs_ot_env(ispin)%ener_gx(:))
1105 0 : tr_xnew_gx = tr_xnew_gx + tmp
1106 : END DO
1107 : END IF
1108 42118 : overlap = (tr_xnew_gx - tr_xold_gx)
1109 : ! OK, bad luck, take a SD step along the preconditioned gradient
1110 42118 : 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 832 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
1118 832 : IF (do_ks) THEN
1119 1774 : DO ispin = 1, nspin
1120 942 : 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 942 : 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 942 : 1.0_dp, 1.0_dp)
1127 1774 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1128 28 : 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 28 : 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 28 : 1.0_dp, 1.0_dp)
1135 : END IF
1136 : END DO
1137 : END IF
1138 832 : IF (do_ener) THEN
1139 0 : DO ispin = 1, nspin
1140 0 : qs_ot_env(ispin)%ener_x(:) = 0._dp
1141 0 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_e(j, :)
1142 0 : 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 42118 : CALL timestop(handle)
1148 :
1149 42118 : 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 178 : 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 178 : 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 178 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f, x
1167 178 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: G, S
1168 : TYPE(cp_logger_type), POINTER :: logger
1169 :
1170 178 : eta = qs_ot_env(1)%settings%broyden_eta
1171 178 : omega = qs_ot_env(1)%settings%broyden_omega
1172 178 : sigma_dec = qs_ot_env(1)%settings%broyden_sigma_decrease
1173 178 : sigma_min = qs_ot_env(1)%settings%broyden_sigma_min
1174 178 : forget_history = qs_ot_env(1)%settings%broyden_forget_history
1175 178 : adaptive_sigma = qs_ot_env(1)%settings%broyden_adaptive_sigma
1176 178 : enable_flip = qs_ot_env(1)%settings%broyden_enable_flip
1177 :
1178 178 : do_ks = qs_ot_env(1)%settings%ks
1179 178 : do_ener = qs_ot_env(1)%settings%do_ener
1180 :
1181 178 : beta = qs_ot_env(1)%settings%broyden_beta
1182 178 : gamma = qs_ot_env(1)%settings%broyden_gamma
1183 178 : IF (adaptive_sigma) THEN
1184 178 : IF (qs_ot_env(1)%broyden_adaptive_sigma .LT. 0.0_dp) THEN
1185 6 : 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 0 : sigma = qs_ot_env(1)%settings%broyden_sigma
1191 : END IF
1192 :
1193 : ! simplify our life....
1194 178 : IF (do_ener .OR. .NOT. do_ks .OR. qs_ot_env(1)%settings%do_rotation) THEN
1195 0 : CPABORT("Not yet implemented")
1196 : END IF
1197 : !
1198 178 : nspin = SIZE(qs_ot_env)
1199 :
1200 178 : diis_m = qs_ot_env(1)%settings%diis_m
1201 :
1202 178 : IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
1203 84 : diis_bound = qs_ot_env(1)%diis_iter + 1
1204 : ELSE
1205 94 : diis_bound = diis_m
1206 : END IF
1207 :
1208 : ! We want x:s, f:s and one random vector
1209 178 : k = 2*diis_bound + 1
1210 712 : ALLOCATE (S(k, k))
1211 712 : ALLOCATE (G(k, k))
1212 534 : ALLOCATE (f(k))
1213 534 : ALLOCATE (x(k))
1214 534 : ALLOCATE (circ_index(diis_bound))
1215 30882 : G = 0.0
1216 2304 : DO i = 1, k
1217 2304 : G(i, i) = sigma
1218 : END DO
1219 30882 : S = 0.0
1220 :
1221 178 : j = MOD(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
1222 :
1223 356 : DO ispin = 1, nspin
1224 356 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
1225 : END DO
1226 :
1227 178 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
1228 178 : qs_ot_env(1)%gnorm = 0.0_dp
1229 356 : DO ispin = 1, nspin
1230 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
1231 178 : 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 178 : tmp)
1234 356 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1235 : END DO
1236 178 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
1237 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1238 : END IF
1239 356 : DO ispin = 1, nspin
1240 356 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -1.0_dp)
1241 : END DO
1242 : ELSE
1243 0 : qs_ot_env(1)%gnorm = 0.0_dp
1244 0 : DO ispin = 1, nspin
1245 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
1246 0 : 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 0 : qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-1.0_dp)
1249 : END DO
1250 : END IF
1251 :
1252 178 : k = 0
1253 : n = 0
1254 178 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
1255 356 : DO ispin = 1, nspin
1256 178 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
1257 356 : k = k + itmp
1258 : END DO
1259 :
1260 : ! Handling the case of no free variables to optimize
1261 178 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) /= 0) THEN
1262 178 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8)))
1263 178 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
1264 : ELSE
1265 0 : qs_ot_env(1)%delta = 0.0_dp
1266 0 : qs_ot_env(1)%gradient = 0.0_dp
1267 : END IF
1268 :
1269 178 : IF (diis_bound == diis_m) THEN
1270 832 : DO i = 1, diis_bound
1271 832 : circ_index(i) = MOD(j + i - 1, diis_m) + 1
1272 : END DO
1273 : ELSE
1274 320 : DO i = 1, diis_bound
1275 320 : circ_index(i) = i
1276 : END DO
1277 : END IF
1278 :
1279 30882 : S = 0.0_dp
1280 356 : DO ispin = 1, nspin
1281 178 : CALL dbcsr_init_random(qs_ot_env(ispin)%matrix_x)
1282 1152 : DO i = 1, diis_bound
1283 : CALL dbcsr_dot( &
1284 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1285 974 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, tmp)
1286 974 : S(i, i) = S(i, i) + tmp
1287 : CALL dbcsr_dot( &
1288 : qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1289 974 : qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, tmp)
1290 974 : 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 974 : qs_ot_env(ispin)%matrix_x, tmp)
1294 974 : S(i, 2*diis_bound + 1) = S(i, 2*diis_bound + 1) + tmp
1295 974 : 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 974 : qs_ot_env(ispin)%matrix_x, tmp)
1299 974 : S(i + diis_bound, 2*diis_bound + 1) = S(i + diis_bound, 2*diis_bound + 1) + tmp
1300 974 : S(i + diis_bound, 2*diis_bound + 1) = S(2*diis_bound + 1, diis_bound + i)
1301 3550 : 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 2576 : tmp)
1306 2576 : S(i, k) = S(i, k) + tmp
1307 2576 : 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 2576 : tmp)
1312 2576 : S(diis_bound + i, diis_bound + k) = S(diis_bound + i, diis_bound + k) + tmp
1313 3550 : S(diis_bound + k, diis_bound + i) = S(diis_bound + i, diis_bound + k)
1314 : END DO
1315 7278 : DO k = 1, diis_bound
1316 : CALL dbcsr_dot( &
1317 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1318 6126 : qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix, tmp)
1319 6126 : S(i, k + diis_bound) = S(i, k + diis_bound) + tmp
1320 7100 : S(k + diis_bound, i) = S(i, k + diis_bound)
1321 : END DO
1322 : END DO
1323 178 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_x, tmp)
1324 356 : 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 178 : k = 2*diis_bound + 1
1329 178 : tmp = SQRT(S(k, k))
1330 2304 : S(k, :) = S(k, :)/tmp
1331 2304 : S(:, k) = S(:, k)/tmp
1332 :
1333 178 : IF (diis_bound .GT. 1) THEN
1334 164 : tmp = 0.0_dp
1335 164 : tmp2 = 0.0_dp
1336 164 : i = diis_bound
1337 328 : 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 164 : tmp)
1343 164 : 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 164 : tmp)
1348 164 : 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 164 : tmp)
1353 164 : 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 164 : tmp)
1358 328 : tmp2 = tmp2 + tmp
1359 : END DO
1360 164 : qs_ot_env(1)%c_broy(i - 1) = tmp2
1361 : END IF
1362 :
1363 178 : qs_ot_env(1)%energy_h(j) = qs_ot_env(1)%etotal
1364 :
1365 : ! If we went uphill, do backtracking line search
1366 1152 : i = MINLOC(qs_ot_env(1)%energy_h(1:diis_bound), dim=1)
1367 178 : IF (i .NE. j) THEN
1368 26 : sigma = sigma_dec*sigma
1369 26 : qs_ot_env(1)%OT_METHOD_FULL = "OT BTRK"
1370 52 : DO ispin = 1, nspin
1371 26 : 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 26 : 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 52 : alpha_scalar=1.0_dp, beta_scalar=gamma)
1378 : END DO
1379 : ELSE
1380 : ! Construct G
1381 792 : DO i = 2, diis_bound
1382 9400 : f = 0.0
1383 9400 : x = 0.0
1384 : ! f is df_i
1385 640 : x(i) = 1.0
1386 640 : x(i - 1) = -1.0
1387 : ! x is dx_i
1388 640 : f(diis_bound + i) = 1.0
1389 640 : f(diis_bound + i - 1) = -1.0
1390 640 : tmp = 1.0_dp
1391 : ! We want a pos def Hessian
1392 640 : IF (enable_flip) THEN
1393 640 : IF (qs_ot_env(1)%c_broy(i - 1) .GT. 0) THEN
1394 : !qs_ot_env(1)%OT_METHOD_FULL="OT FLIP"
1395 2 : tmp = -1.0_dp
1396 : END IF
1397 : END IF
1398 :
1399 : ! get dx-Gdf
1400 400200 : x(:) = tmp*x - MATMUL(G, f)
1401 : ! dfSdf
1402 : ! we calculate matmul(S, f) twice. They're small...
1403 400200 : tmp = DOT_PRODUCT(f, MATMUL(S, f))
1404 : ! NOTE THAT S IS SYMMETRIC !!!
1405 400200 : f(:) = MATMUL(S, f)/tmp
1406 : ! the spread is an outer vector product
1407 133552 : G(:, :) = G + SPREAD(x, dim=2, ncopies=SIZE(f))*SPREAD(f, dim=1, ncopies=SIZE(x))
1408 : END DO
1409 1888 : f = 0.0_dp
1410 152 : f(2*diis_bound) = 1.0_dp
1411 96272 : x(:) = -beta*MATMUL(G, f)
1412 :
1413 : ! OK, add the vectors now, this sums up to the proposed step
1414 304 : DO ispin = 1, nspin
1415 152 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1416 944 : 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 944 : alpha_scalar=1.0_dp, beta_scalar=-x(i + diis_bound))
1420 : END DO
1421 1096 : 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 944 : alpha_scalar=1.0_dp, beta_scalar=x(i))
1425 : END DO
1426 : END DO
1427 :
1428 152 : IF (adaptive_sigma) THEN
1429 152 : tmp = new_sigma(G, S, diis_bound)
1430 : !tmp = tmp * qs_ot_env ( 1 ) % settings % broyden_sigma
1431 152 : tmp = tmp*eta
1432 152 : sigma = MIN(omega*sigma, tmp)
1433 : END IF
1434 :
1435 : ! compute the inner product of direction of the step and gradient
1436 152 : tmp = 0.0_dp
1437 304 : DO ispin = 1, nspin
1438 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
1439 : qs_ot_env(ispin)%matrix_x, &
1440 152 : tmp2)
1441 304 : tmp = tmp + tmp2
1442 : END DO
1443 :
1444 304 : DO ispin = 1, nspin
1445 : ! if the direction of the step is not in direction of the gradient,
1446 : ! change step sign
1447 304 : IF (tmp .GE. 0.0_dp) THEN
1448 8 : qs_ot_env(1)%OT_METHOD_FULL = "OT TURN"
1449 8 : IF (forget_history) THEN
1450 0 : qs_ot_env(1)%diis_iter = 0
1451 : END IF
1452 8 : 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 8 : 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 144 : 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 178 : DEALLOCATE (S, G, f, x, circ_index)
1466 :
1467 : ! update for next round
1468 178 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1469 178 : qs_ot_env(1)%broyden_adaptive_sigma = MAX(sigma, sigma_min)
1470 :
1471 356 : END SUBROUTINE ot_broyden_step
1472 :
1473 : ! **************************************************************************************************
1474 : !> \brief ...
1475 : !> \param G ...
1476 : !> \param S ...
1477 : !> \param n ...
1478 : !> \return ...
1479 : ! **************************************************************************************************
1480 152 : 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 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H
1492 :
1493 608 : ALLOCATE (H(n, n))
1494 152 : CALL hess_G(G, S, H, n)
1495 456 : ALLOCATE (eigv(n))
1496 152 : 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 1864 : 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 152 : DEALLOCATE (H, eigv)
1511 152 : END FUNCTION new_sigma
1512 :
1513 : ! **************************************************************************************************
1514 : !> \brief ...
1515 : !> \param G ...
1516 : !> \param S ...
1517 : !> \param H ...
1518 : !> \param n ...
1519 : ! **************************************************************************************************
1520 152 : 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 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: v
1533 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Q
1534 :
1535 152 : i = SIZE(G, 1)
1536 152 : k = SIZE(H, 1)
1537 608 : ALLOCATE (Q(i, k))
1538 456 : ALLOCATE (v(i))
1539 5796 : H = 0.0_dp
1540 11440 : Q = 0.0_dp
1541 :
1542 1888 : Q(:, 1) = 1.0_dp
1543 28392 : tmp = SQRT(DOT_PRODUCT(Q(:, 1), MATMUL(S, Q(:, 1))))
1544 11440 : Q(:, :) = Q(:, :)/tmp
1545 :
1546 926 : DO i = 1, k
1547 165440 : v(:) = MATMUL(G, Q(:, i))
1548 3530 : DO j = 1, i
1549 669822 : H(j, i) = DOT_PRODUCT(Q(:, j), MATMUL(S, v))
1550 41756 : v(:) = v - H(j, i)*Q(:, j)
1551 : END DO
1552 926 : IF (i .LT. k) THEN
1553 150562 : tmp = DOT_PRODUCT(v, MATMUL(S, v))
1554 634 : IF (tmp .LE. 0.0_dp) THEN
1555 4 : n = i
1556 4 : EXIT
1557 : END IF
1558 630 : tmp = SQRT(tmp)
1559 : ! Lucky breakdown
1560 630 : IF (ABS(tmp) .LT. 1e-9_dp) THEN
1561 2 : n = i
1562 2 : EXIT
1563 : END IF
1564 628 : H(i + 1, i) = tmp
1565 9208 : Q(:, i + 1) = v/H(i + 1, i)
1566 : END IF
1567 : END DO
1568 :
1569 152 : DEALLOCATE (Q, v)
1570 152 : END SUBROUTINE hess_G
1571 :
1572 5456 : END MODULE qs_ot_minimizer
|
