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