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