Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
44 : CONTAINS
45 : !
46 : ! the minimizer interface
47 : ! should present all possible modes of minimization
48 : ! these include CG SD DIIS
49 : !
50 : !
51 : ! IN the case of nspin != 1 we have a gradient that is distributed over different qs_ot_env.
52 : ! still things remain basically the same, since there are no constraints between the different qs_ot_env
53 : ! we only should take care that the various scalar products are taken over the full vectors.
54 : ! all the information needed and collected can be stored in the fist qs_ot_env only
55 : ! (indicating that the data type for the gradient/position and minization should be separated)
56 : !
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param qs_ot_env ...
60 : !> \param matrix_hc ...
61 : ! **************************************************************************************************
62 82346 : 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 82346 : CALL timeset(routineN, handle)
73 :
74 82346 : nspin = SIZE(qs_ot_env)
75 :
76 82346 : do_ks = qs_ot_env(1)%settings%ks
77 82346 : do_ener = qs_ot_env(1)%settings%do_ener
78 :
79 82346 : qs_ot_env(1)%OT_METHOD_FULL = ""
80 :
81 : ! compute the gradient for the variables x
82 82346 : IF (.NOT. qs_ot_env(1)%energy_only) THEN
83 62356 : qs_ot_env(1)%gradient = 0.0_dp
84 132201 : DO ispin = 1, nspin
85 69845 : IF (do_ks) THEN
86 136694 : 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 66849 : 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 2996 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin))
95 : CASE DEFAULT
96 69845 : CPABORT("ALGORITHM NYI")
97 : END SELECT
98 : END IF
99 : ! and also the gradient along the direction
100 132201 : IF (qs_ot_env(1)%use_dx) THEN
101 65101 : IF (do_ks) THEN
102 65101 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
103 65101 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
104 65101 : IF (qs_ot_env(1)%settings%do_rotation) THEN
105 1630 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
106 1630 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + 0.5_dp*tmp
107 : END IF
108 : END IF
109 65101 : IF (do_ener) THEN
110 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
111 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient + tmp
112 : END IF
113 : ELSE
114 4744 : IF (do_ks) THEN
115 4744 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
116 4744 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
117 4744 : IF (qs_ot_env(1)%settings%do_rotation) THEN
118 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
119 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - 0.5_dp*tmp
120 : END IF
121 : END IF
122 4744 : IF (do_ener) THEN
123 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
124 0 : qs_ot_env(1)%gradient = qs_ot_env(1)%gradient - tmp
125 : END IF
126 : END IF
127 : END DO
128 : END IF
129 :
130 123940 : SELECT CASE (qs_ot_env(1)%settings%OT_METHOD)
131 : CASE ("CG")
132 41594 : IF (current_point_is_fine(qs_ot_env)) THEN
133 21698 : qs_ot_env(1)%OT_METHOD_FULL = "OT CG"
134 21698 : CALL ot_new_cg_direction(qs_ot_env)
135 21698 : qs_ot_env(1)%line_search_count = 0
136 : ELSE
137 19896 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
138 : END IF
139 41594 : CALL do_line_search(qs_ot_env)
140 : CASE ("SD")
141 132 : IF (current_point_is_fine(qs_ot_env)) THEN
142 38 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
143 38 : CALL ot_new_sd_direction(qs_ot_env)
144 38 : qs_ot_env(1)%line_search_count = 0
145 : ELSE
146 94 : qs_ot_env(1)%OT_METHOD_FULL = "OT LS"
147 : END IF
148 132 : CALL do_line_search(qs_ot_env)
149 : CASE ("DIIS")
150 40444 : qs_ot_env(1)%OT_METHOD_FULL = "OT DIIS"
151 40444 : CALL ot_diis_step(qs_ot_env)
152 : CASE ("BROY")
153 176 : qs_ot_env(1)%OT_METHOD_FULL = "OT BROY"
154 176 : CALL ot_broyden_step(qs_ot_env)
155 : CASE DEFAULT
156 82346 : CPABORT("OT_METHOD NYI")
157 : END SELECT
158 :
159 82346 : CALL timestop(handle)
160 :
161 82346 : 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 41726 : 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 41726 : res = .FALSE.
177 :
178 : ! only if we have a gradient it can be fine
179 41726 : IF (.NOT. qs_ot_env(1)%energy_only) THEN
180 :
181 : ! we have not yet started with the line search
182 21736 : IF (qs_ot_env(1)%line_search_count .EQ. 0) THEN
183 41726 : res = .TRUE.
184 : RETURN
185 : END IF
186 :
187 18798 : IF (qs_ot_env(1)%line_search_might_be_done) THEN
188 : ! here we put the more complicated logic later
189 18798 : res = .TRUE.
190 18798 : 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 41726 : SUBROUTINE do_line_search(qs_ot_env)
205 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
206 :
207 42154 : SELECT CASE (qs_ot_env(1)%settings%line_search_method)
208 : CASE ("GOLD")
209 428 : CALL do_line_search_gold(qs_ot_env)
210 : CASE ("3PNT")
211 1664 : CALL do_line_search_3pnt(qs_ot_env)
212 : CASE ("2PNT")
213 39394 : CALL do_line_search_2pnt(qs_ot_env)
214 : CASE ("ADPT")
215 230 : CALL do_line_search_adapt(qs_ot_env)
216 : CASE ("NONE")
217 10 : CALL do_line_search_none(qs_ot_env)
218 : CASE DEFAULT
219 41726 : CPABORT("NYI")
220 : END SELECT
221 41726 : 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 41726 : 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 41726 : CALL timeset(routineN, handle)
240 :
241 41726 : nspin = SIZE(qs_ot_env)
242 :
243 41726 : do_ks = qs_ot_env(1)%settings%ks
244 41726 : 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 41726 : IF (qs_ot_env(1)%use_dx) THEN
249 41638 : IF (do_ks) THEN
250 91471 : DO ispin = 1, nspin
251 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_dx, &
252 49833 : alpha_scalar=1.0_dp, beta_scalar=ds)
253 91471 : 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 2644 : alpha_scalar=1.0_dp, beta_scalar=ds)
256 : END IF
257 : END DO
258 : END IF
259 41638 : IF (do_ener) THEN
260 0 : DO ispin = 1, nspin
261 0 : 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 88 : IF (do_ks) THEN
266 176 : DO ispin = 1, nspin
267 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_gx, &
268 88 : alpha_scalar=1.0_dp, beta_scalar=-ds)
269 176 : 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 0 : alpha_scalar=1.0_dp, beta_scalar=-ds)
272 : END IF
273 : END DO
274 : END IF
275 88 : IF (do_ener) THEN
276 0 : DO ispin = 1, nspin
277 0 : 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 41726 : CALL timestop(handle)
282 41726 : 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 428 : 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 428 : CALL timeset(routineN, handle)
300 :
301 428 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
302 428 : count = qs_ot_env(1)%line_search_count
303 428 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
304 428 : qs_ot_env(1)%energy_only = .TRUE.
305 :
306 428 : 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 0 : CPABORT("MAX ITER EXCEEDED : FATAL")
310 : END IF
311 :
312 428 : IF (qs_ot_env(1)%line_search_count .EQ. 1) THEN
313 54 : qs_ot_env(1)%line_search_left = 1
314 54 : qs_ot_env(1)%line_search_right = 0
315 54 : qs_ot_env(1)%line_search_mid = 1
316 54 : qs_ot_env(1)%ot_pos(1) = 0.0_dp
317 54 : qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
318 54 : qs_ot_env(1)%ot_pos(2) = qs_ot_env(1)%ds_min/gold_sec
319 : ELSE
320 374 : 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 374 : IF (qs_ot_env(1)%line_search_right .EQ. 0) THEN ! we do not yet have the right bracket
325 64 : IF (qs_ot_env(1)%ot_energy(count - 1) .LT. qs_ot_env(1)%ot_energy(count)) THEN
326 42 : 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 42 : qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid))*gold_sec
330 : ELSE
331 22 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
332 22 : qs_ot_env(1)%line_search_mid = count
333 22 : 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 310 : IF (qs_ot_env(1)%ot_pos(count) .LT. qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
338 126 : IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
339 42 : qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
340 42 : qs_ot_env(1)%line_search_mid = count
341 : ELSE
342 84 : qs_ot_env(1)%line_search_left = count
343 : END IF
344 : ELSE
345 184 : IF (qs_ot_env(1)%ot_energy(count) .LT. qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
346 76 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
347 76 : qs_ot_env(1)%line_search_mid = count
348 : ELSE
349 108 : 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 310 : 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 162 : - 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 148 : - 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 310 : 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 42 : qs_ot_env(1)%energy_only = .FALSE.
375 42 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
376 : END IF
377 : END IF
378 : END IF
379 428 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
380 428 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
381 :
382 428 : CALL take_step(ds, qs_ot_env)
383 :
384 428 : CALL timestop(handle)
385 :
386 428 : END SUBROUTINE do_line_search_gold
387 :
388 : ! **************************************************************************************************
389 : !> \brief ...
390 : !> \param qs_ot_env ...
391 : ! **************************************************************************************************
392 1664 : 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 1664 : CALL timeset(routineN, handle)
403 :
404 1664 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
405 1664 : qs_ot_env(1)%energy_only = .TRUE.
406 :
407 : ! a three point interpolation based on the energy
408 1664 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
409 1664 : count = qs_ot_env(1)%line_search_count
410 1664 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
411 580 : SELECT CASE (count)
412 : CASE (1)
413 580 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
414 580 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*0.8_dp
415 : CASE (2)
416 542 : IF (qs_ot_env(1)%OT_energy(count) .GT. qs_ot_env(1)%OT_energy(count - 1)) THEN
417 70 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*0.5_dp
418 : ELSE
419 472 : qs_ot_env(1)%OT_pos(count + 1) = qs_ot_env(1)%ds_min*1.4_dp
420 : END IF
421 : CASE (3)
422 542 : xa = qs_ot_env(1)%OT_pos(1)
423 542 : xb = qs_ot_env(1)%OT_pos(2)
424 542 : xc = qs_ot_env(1)%OT_pos(3)
425 542 : fa = qs_ot_env(1)%OT_energy(1)
426 542 : fb = qs_ot_env(1)%OT_energy(2)
427 542 : fc = qs_ot_env(1)%OT_energy(3)
428 542 : nom = (xb - xa)**2*(fb - fc) - (xb - xc)**2*(fb - fa)
429 542 : denom = (xb - xa)*(fb - fc) - (xb - xc)*(fb - fa)
430 542 : IF (ABS(denom) .LE. 1.0E-18_dp*MAX(ABS(fb - fc), ABS(fb - fa))) THEN
431 : pos = xb
432 : ELSE
433 542 : 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 542 : (pos - xc)*(pos - xa)*fb/((xb - xc)*(xb - xa))
438 542 : 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 2710 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:3))*4.0_dp))
442 : ELSE ! just take an extended step
443 0 : qs_ot_env(1)%OT_pos(count + 1) = MAXVAL(qs_ot_env(1)%OT_pos(1:3))*2.0_dp
444 : END IF
445 542 : qs_ot_env(1)%energy_only = .FALSE.
446 542 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
447 : CASE DEFAULT
448 1664 : CPABORT("NYI")
449 : END SELECT
450 1664 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
451 1664 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
452 :
453 1664 : CALL take_step(ds, qs_ot_env)
454 :
455 1664 : CALL timestop(handle)
456 :
457 1664 : END SUBROUTINE do_line_search_3pnt
458 :
459 : ! **************************************************************************************************
460 : !> \brief ...
461 : !> \param qs_ot_env ...
462 : ! **************************************************************************************************
463 39394 : 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 39394 : CALL timeset(routineN, handle)
473 :
474 39394 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
475 39394 : qs_ot_env(1)%energy_only = .TRUE.
476 :
477 : ! a three point interpolation based on the energy
478 39394 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
479 39394 : count = qs_ot_env(1)%line_search_count
480 39394 : qs_ot_env(1)%ot_energy(count) = qs_ot_env(1)%etotal
481 21024 : SELECT CASE (count)
482 : CASE (1)
483 21024 : qs_ot_env(1)%ot_pos(count) = 0.0_dp
484 21024 : qs_ot_env(1)%ot_grad(count) = qs_ot_env(1)%gradient
485 21024 : qs_ot_env(1)%ot_pos(count + 1) = qs_ot_env(1)%ds_min*1.0_dp
486 : CASE (2)
487 18370 : x0 = 0.0_dp
488 18370 : c = qs_ot_env(1)%ot_energy(1)
489 18370 : b = qs_ot_env(1)%ot_grad(1)
490 18370 : x1 = qs_ot_env(1)%ot_pos(2)
491 18370 : a = (qs_ot_env(1)%ot_energy(2) - b*x1 - c)/(x1**2)
492 18370 : IF (a .LE. 0.0_dp) a = 1.0E-15_dp
493 18370 : pos = -b/(2.0_dp*a)
494 18370 : val = a*pos**2 + b*pos + c
495 18370 : qs_ot_env(1)%energy_only = .FALSE.
496 18370 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
497 18370 : 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 73232 : MIN(pos, MAXVAL(qs_ot_env(1)%OT_pos(1:2))*4.0_dp))
502 : ELSE ! just take an extended step
503 248 : 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 39394 : CPABORT("NYI")
507 : END SELECT
508 39394 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
509 39394 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
510 :
511 39394 : CALL take_step(ds, qs_ot_env)
512 :
513 39394 : CALL timestop(handle)
514 :
515 39394 : END SUBROUTINE do_line_search_2pnt
516 :
517 : ! **************************************************************************************************
518 : !> \brief ...
519 : !> \param qs_ot_env ...
520 : ! **************************************************************************************************
521 230 : 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 230 : CALL timeset(routineN, handle)
534 :
535 230 : qs_ot_env(1)%line_search_count = qs_ot_env(1)%line_search_count + 1
536 230 : count = qs_ot_env(1)%line_search_count
537 230 : qs_ot_env(1)%line_search_might_be_done = .FALSE.
538 230 : qs_ot_env(1)%energy_only = .TRUE.
539 :
540 230 : 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 0 : CPABORT("MAX ITER EXCEEDED : FATAL")
544 : END IF
545 :
546 : ! Perform an adaptive linesearch
547 230 : IF (qs_ot_env(1)%line_search_count == 1) THEN
548 68 : qs_ot_env(1)%line_search_left = 1
549 68 : qs_ot_env(1)%line_search_right = 0
550 68 : qs_ot_env(1)%line_search_mid = 1
551 68 : qs_ot_env(1)%ot_Pos(1) = 0.0_dp
552 68 : qs_ot_env(1)%ot_energy(1) = qs_ot_env(1)%etotal
553 68 : qs_ot_env(1)%ot_Pos(2) = qs_ot_env(1)%ds_min*grow_factor
554 : ELSE
555 162 : 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 162 : IF (qs_ot_env(1)%line_search_right == 0) THEN ! we do not yet have the right bracket
560 88 : IF (qs_ot_env(1)%ot_energy(count - 1) < qs_ot_env(1)%ot_energy(count)) THEN
561 58 : 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 58 : qs_ot_env(1)%ot_Pos(qs_ot_env(1)%line_search_mid))*shrink_factor
565 : ELSE
566 : ! expand further
567 30 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
568 30 : qs_ot_env(1)%line_search_mid = count
569 30 : 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 74 : IF (qs_ot_env(1)%ot_pos(count) < qs_ot_env(1)%ot_pos(qs_ot_env(1)%line_search_mid)) THEN
574 0 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
575 0 : qs_ot_env(1)%line_search_right = qs_ot_env(1)%line_search_mid
576 0 : qs_ot_env(1)%line_search_mid = count
577 : ELSE
578 0 : qs_ot_env(1)%line_search_left = count
579 : END IF
580 : ELSE
581 74 : IF (qs_ot_env(1)%ot_energy(count) < qs_ot_env(1)%ot_energy(qs_ot_env(1)%line_search_mid)) THEN
582 36 : qs_ot_env(1)%line_search_left = qs_ot_env(1)%line_search_mid
583 36 : qs_ot_env(1)%line_search_mid = count
584 : ELSE
585 38 : qs_ot_env(1)%line_search_right = count
586 : END IF
587 : END IF
588 74 : il = qs_ot_env(1)%line_search_left
589 74 : im = qs_ot_env(1)%line_search_mid
590 74 : ir = qs_ot_env(1)%line_search_right
591 74 : xl = qs_ot_env(1)%OT_pos(il)
592 74 : xm = qs_ot_env(1)%OT_pos(im)
593 74 : xr = qs_ot_env(1)%OT_pos(ir)
594 74 : el = qs_ot_env(1)%ot_energy(il)
595 74 : em = qs_ot_env(1)%ot_energy(im)
596 74 : er = qs_ot_env(1)%ot_energy(ir)
597 74 : IF (em < el) THEN
598 58 : IF (er < em) THEN
599 : !extend search
600 0 : qs_ot_env(1)%ot_Pos(count + 1) = qs_ot_env(1)%ot_Pos(ir)*grow_factor
601 : ELSE
602 : ! Cramer's rule
603 58 : denom = (xl - xm)*(xl - xr)*(xm - xr)
604 58 : a = (xr*(em - el) + xm*(el - er) + xl*(er - em))/denom
605 58 : b = (xr**2*(el - em) + xm**2*(er - el) + xl**2*(em - er))/denom
606 58 : c = (xm*xr*(xm - xr)*el + xr*xl*(xr - xl)*em + xr*xm*(xr - xm)*er)/denom
607 :
608 58 : IF (ABS(a) /= 0.0_dp) THEN
609 58 : step_size = -b/(2.0_dp*a)
610 : ELSE
611 : step_size = 0.0_dp
612 : END IF
613 58 : CPASSERT(step_size >= 0.0_dp)
614 58 : qs_ot_env(1)%ot_Pos(count + 1) = step_size
615 58 : qs_ot_env(1)%line_search_might_be_done = .TRUE.
616 58 : 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 16 : (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 230 : ds = qs_ot_env(1)%OT_pos(count + 1) - qs_ot_env(1)%OT_pos(count)
627 230 : qs_ot_env(1)%ds_min = qs_ot_env(1)%OT_pos(count + 1)
628 :
629 230 : CALL take_step(ds, qs_ot_env)
630 :
631 230 : CALL timestop(handle)
632 :
633 230 : END SUBROUTINE do_line_search_adapt
634 :
635 : ! **************************************************************************************************
636 : !> \brief ...
637 : !> \param qs_ot_env ...
638 : ! **************************************************************************************************
639 10 : SUBROUTINE do_line_search_none(qs_ot_env)
640 : TYPE(qs_ot_type), DIMENSION(:), POINTER :: qs_ot_env
641 :
642 10 : CALL take_step(qs_ot_env(1)%ds_min, qs_ot_env)
643 :
644 10 : 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 38 : 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 38 : CALL timeset(routineN, handle)
666 :
667 : !***SCP
668 :
669 38 : nspin = SIZE(qs_ot_env)
670 38 : logger => cp_get_default_logger()
671 38 : do_ks = qs_ot_env(1)%settings%ks
672 38 : do_ener = qs_ot_env(1)%settings%do_ener
673 :
674 38 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
675 26 : IF (.NOT. qs_ot_env(1)%use_dx) CPABORT("use dx")
676 26 : qs_ot_env(1)%gnorm = 0.0_dp
677 26 : IF (do_ks) THEN
678 52 : DO ispin = 1, nspin
679 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
680 26 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx)
681 26 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
682 52 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
683 : END DO
684 26 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
685 0 : logger => cp_get_default_logger()
686 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
687 : END IF
688 52 : DO ispin = 1, nspin
689 52 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_dx, -1.0_dp)
690 : END DO
691 26 : IF (qs_ot_env(1)%settings%do_rotation) THEN
692 0 : DO ispin = 1, nspin
693 : ! right now no preconditioner yet
694 0 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_dx, qs_ot_env(ispin)%rot_mat_gx)
695 0 : 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 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
698 : END DO
699 0 : DO ispin = 1, nspin
700 0 : CALL dbcsr_scale(qs_ot_env(ispin)%rot_mat_dx, -1.0_dp)
701 : END DO
702 : END IF
703 : END IF
704 26 : IF (do_ener) THEN
705 0 : DO ispin = 1, nspin
706 0 : qs_ot_env(ispin)%ener_dx = qs_ot_env(ispin)%ener_gx
707 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_dx, qs_ot_env(ispin)%ener_gx)
708 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
709 0 : qs_ot_env(ispin)%ener_dx = -qs_ot_env(ispin)%ener_dx
710 : END DO
711 : END IF
712 : ELSE
713 12 : qs_ot_env(1)%gnorm = 0.0_dp
714 12 : IF (do_ks) THEN
715 24 : DO ispin = 1, nspin
716 12 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
717 24 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
718 : END DO
719 12 : IF (qs_ot_env(1)%settings%do_rotation) THEN
720 0 : DO ispin = 1, nspin
721 0 : 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 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
724 : END DO
725 : END IF
726 : END IF
727 12 : IF (do_ener) THEN
728 0 : DO ispin = 1, nspin
729 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
730 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
731 : END DO
732 : END IF
733 : END IF
734 :
735 38 : k = 0
736 38 : n = 0
737 38 : nener = 0
738 38 : IF (do_ks) THEN
739 38 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
740 76 : DO ispin = 1, nspin
741 38 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
742 76 : k = k + itmp
743 : END DO
744 : END IF
745 38 : IF (do_ener) THEN
746 0 : DO ispin = 1, nspin
747 0 : 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 38 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
752 36 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
753 36 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
754 : ELSE
755 2 : qs_ot_env(1)%delta = 0.0_dp
756 2 : qs_ot_env(1)%gradient = 0.0_dp
757 : END IF
758 :
759 38 : CALL timestop(handle)
760 :
761 38 : 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 21698 : 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 21698 : CALL timeset(routineN, handle)
783 :
784 21698 : nspin = SIZE(qs_ot_env)
785 21698 : logger => cp_get_default_logger()
786 :
787 21698 : do_ks = qs_ot_env(1)%settings%ks
788 21698 : do_ener = qs_ot_env(1)%settings%do_ener
789 :
790 21698 : gnorm_cross = 0.0_dp
791 21698 : IF (do_ks) THEN
792 47895 : DO ispin = 1, nspin
793 26197 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
794 47895 : gnorm_cross = gnorm_cross + tmp
795 : END DO
796 21698 : IF (qs_ot_env(1)%settings%do_rotation) THEN
797 2458 : DO ispin = 1, nspin
798 1264 : 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 2458 : gnorm_cross = gnorm_cross + 0.5_dp*tmp
801 : END DO
802 : END IF
803 : END IF
804 21698 : IF (do_ener) THEN
805 0 : DO ispin = 1, nspin
806 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
807 0 : gnorm_cross = gnorm_cross + tmp
808 : END DO
809 : END IF
810 :
811 21698 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
812 :
813 36219 : DO ispin = 1, nspin
814 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
815 36219 : qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
816 : END DO
817 16100 : qs_ot_env(1)%gnorm = 0.0_dp
818 16100 : IF (do_ks) THEN
819 36219 : DO ispin = 1, nspin
820 20119 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old, tmp)
821 36219 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
822 : END DO
823 16100 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
824 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
825 : END IF
826 36219 : DO ispin = 1, nspin
827 36219 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx_old)
828 : END DO
829 16100 : IF (qs_ot_env(1)%settings%do_rotation) THEN
830 2458 : DO ispin = 1, nspin
831 : ! right now no preconditioner yet
832 1264 : CALL dbcsr_copy(qs_ot_env(ispin)%rot_mat_gx_old, qs_ot_env(ispin)%rot_mat_gx)
833 1264 : 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 2458 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
836 : END DO
837 2458 : DO ispin = 1, nspin
838 2458 : 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 16100 : IF (do_ener) THEN
843 0 : DO ispin = 1, nspin
844 0 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
845 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx_old)
846 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
847 0 : qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_gx_old
848 : END DO
849 : END IF
850 : ELSE
851 5598 : IF (do_ks) THEN
852 5598 : qs_ot_env(1)%gnorm = 0.0_dp
853 11676 : DO ispin = 1, nspin
854 6078 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
855 6078 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
856 11676 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_gx_old, qs_ot_env(ispin)%matrix_gx)
857 : END DO
858 5598 : IF (qs_ot_env(1)%settings%do_rotation) THEN
859 0 : DO ispin = 1, nspin
860 0 : 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 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
863 0 : 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 5598 : IF (do_ener) THEN
868 0 : DO ispin = 1, nspin
869 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_gx)
870 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
871 0 : qs_ot_env(ispin)%ener_gx_old = qs_ot_env(ispin)%ener_gx
872 : END DO
873 : END IF
874 : END IF
875 :
876 21698 : k = 0
877 21698 : n = 0
878 21698 : nener = 0
879 21698 : IF (do_ks) THEN
880 21698 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
881 47895 : DO ispin = 1, nspin
882 26197 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
883 47895 : k = k + itmp
884 : END DO
885 : END IF
886 21698 : IF (do_ener) THEN
887 0 : DO ispin = 1, nspin
888 0 : 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 21698 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
893 21614 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
894 21614 : beta_pr = (qs_ot_env(1)%gnorm - gnorm_cross)/qs_ot_env(1)%gnorm_old
895 : ELSE
896 84 : qs_ot_env(1)%delta = 0.0_dp
897 84 : beta_pr = 0.0_dp
898 : END IF
899 21698 : beta_pr = MAX(beta_pr, 0.0_dp) ! reset to SD
900 :
901 21698 : test_down = 0.0_dp
902 21698 : IF (do_ks) THEN
903 47895 : DO ispin = 1, nspin
904 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
905 26197 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
906 26197 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_dx, tmp)
907 26197 : test_down = test_down + tmp
908 47895 : 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 1264 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
911 1264 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_dx, tmp)
912 1264 : test_down = test_down + 0.5_dp*tmp
913 : END IF
914 : END DO
915 : END IF
916 21698 : IF (do_ener) THEN
917 0 : DO ispin = 1, nspin
918 0 : qs_ot_env(ispin)%ener_dx = beta_pr*qs_ot_env(ispin)%ener_dx - qs_ot_env(ispin)%ener_gx
919 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx, qs_ot_env(ispin)%ener_dx)
920 0 : test_down = test_down + tmp
921 : END DO
922 : END IF
923 :
924 21698 : IF (test_down .GE. 0.0_dp) THEN ! reset to SD
925 564 : beta_pr = 0.0_dp
926 564 : IF (do_ks) THEN
927 1292 : DO ispin = 1, nspin
928 : CALL dbcsr_add(qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin)%matrix_gx, &
929 728 : alpha_scalar=beta_pr, beta_scalar=-1.0_dp)
930 1292 : IF (qs_ot_env(1)%settings%do_rotation) THEN
931 : CALL dbcsr_add(qs_ot_env(ispin)%rot_mat_dx, &
932 10 : 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 564 : IF (do_ener) THEN
937 0 : DO ispin = 1, nspin
938 0 : 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 21698 : qs_ot_env(1)%gradient = beta_pr*qs_ot_env(1)%gradient - qs_ot_env(1)%gnorm
944 21698 : qs_ot_env(1)%gnorm_old = qs_ot_env(1)%gnorm
945 :
946 21698 : CALL timestop(handle)
947 :
948 21698 : END SUBROUTINE ot_new_cg_direction
949 :
950 : ! **************************************************************************************************
951 : !> \brief ...
952 : !> \param qs_ot_env ...
953 : ! **************************************************************************************************
954 40444 : 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 40444 : CALL timeset(routineN, handle)
966 :
967 40444 : logger => cp_get_default_logger()
968 :
969 40444 : do_ks = qs_ot_env(1)%settings%ks
970 40444 : do_ener = qs_ot_env(1)%settings%do_ener
971 40444 : nspin = SIZE(qs_ot_env)
972 :
973 40444 : diis_m = qs_ot_env(1)%settings%diis_m
974 :
975 40444 : IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
976 25402 : diis_bound = qs_ot_env(1)%diis_iter + 1
977 : ELSE
978 : diis_bound = diis_m
979 : END IF
980 :
981 40444 : 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 40444 : IF (do_ks) THEN
986 83878 : DO ispin = 1, nspin
987 43434 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
988 83878 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
989 366 : 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 40444 : IF (do_ener) THEN
994 0 : DO ispin = 1, nspin
995 0 : qs_ot_env(ispin)%ener_h_x(j, :) = qs_ot_env(ispin)%ener_x(:)
996 : END DO
997 : END IF
998 40444 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
999 35712 : qs_ot_env(1)%gnorm = 0.0_dp
1000 35712 : IF (do_ks) THEN
1001 74414 : DO ispin = 1, nspin
1002 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
1003 38702 : 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 38702 : tmp)
1006 74414 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1007 : END DO
1008 35712 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
1009 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1010 : END IF
1011 74414 : DO ispin = 1, nspin
1012 74414 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -qs_ot_env(1)%ds_min)
1013 : END DO
1014 35712 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1015 732 : DO ispin = 1, nspin
1016 366 : 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 366 : tmp)
1019 732 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + 0.5_dp*tmp
1020 : END DO
1021 732 : DO ispin = 1, nspin
1022 732 : 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 35712 : IF (do_ener) THEN
1027 0 : DO ispin = 1, nspin
1028 0 : qs_ot_env(ispin)%ener_h_e(j, :) = qs_ot_env(ispin)%ener_gx(:)
1029 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_e(j, :), qs_ot_env(ispin)%ener_gx(:))
1030 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1031 0 : 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 4732 : qs_ot_env(1)%gnorm = 0.0_dp
1036 4732 : IF (do_ks) THEN
1037 9464 : DO ispin = 1, nspin
1038 4732 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
1039 4732 : 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 9464 : qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-qs_ot_env(1)%ds_min)
1042 : END DO
1043 4732 : IF (qs_ot_env(1)%settings%do_rotation) THEN
1044 0 : DO ispin = 1, nspin
1045 0 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_gx, qs_ot_env(ispin)%rot_mat_gx, tmp)
1046 0 : 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 0 : 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 4732 : IF (do_ener) THEN
1053 0 : DO ispin = 1, nspin
1054 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_gx(:))
1055 0 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1056 0 : 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 40444 : k = 0
1061 40444 : n = 0
1062 40444 : nener = 0
1063 40444 : IF (do_ks) THEN
1064 40444 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
1065 83878 : DO ispin = 1, nspin
1066 43434 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
1067 83878 : k = k + itmp
1068 : END DO
1069 : END IF
1070 40444 : IF (do_ener) THEN
1071 0 : DO ispin = 1, nspin
1072 0 : 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 40444 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener /= 0) THEN
1077 40418 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8) + nener))
1078 40418 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
1079 : ELSE
1080 26 : qs_ot_env(1)%delta = 0.0_dp
1081 26 : qs_ot_env(1)%gradient = 0.0_dp
1082 : END IF
1083 :
1084 : ! make the diis matrix and solve it
1085 244061 : 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 203617 : qs_ot_env(1)%ls_diis(i, j) = 0.0_dp
1114 203617 : IF (do_ks) THEN
1115 423504 : 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 219887 : tmp)
1119 219887 : qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*tmp
1120 423504 : 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 1976 : tmp)
1124 1976 : qs_ot_env(1)%ls_diis(i, j) = qs_ot_env(1)%ls_diis(i, j) - qs_ot_env(1)%ds_min*0.5_dp*tmp
1125 : END IF
1126 : END DO
1127 : END IF
1128 203617 : IF (do_ener) THEN
1129 0 : DO ispin = 1, nspin
1130 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_gx(:), qs_ot_env(ispin)%ener_h_e(i, :))
1131 0 : 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 203617 : qs_ot_env(1)%ls_diis(j, i) = qs_ot_env(1)%ls_diis(i, j)
1136 203617 : qs_ot_env(1)%ls_diis(i, diis_bound + 1) = 1.0_dp
1137 203617 : qs_ot_env(1)%ls_diis(diis_bound + 1, i) = 1.0_dp
1138 244061 : qs_ot_env(1)%c_diis(i) = 0.0_dp
1139 : END DO
1140 40444 : qs_ot_env(1)%ls_diis(diis_bound + 1, diis_bound + 1) = 0.0_dp
1141 40444 : qs_ot_env(1)%c_diis(diis_bound + 1) = 1.0_dp
1142 : ! put in buffer, dgesv destroys
1143 3075912 : 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 40444 : qs_ot_env(1)%c_diis, diis_m + 1, info)
1147 :
1148 40444 : IF (info .NE. 0) THEN
1149 2 : do_ot_sd = .TRUE.
1150 2 : WRITE (cp_logger_get_default_unit_nr(logger), *) "Singular DIIS matrix"
1151 : ELSE
1152 40442 : do_ot_sd = .FALSE.
1153 40442 : IF (do_ks) THEN
1154 83874 : DO ispin = 1, nspin
1155 : ! OK, add the vectors now
1156 43432 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1157 263305 : 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 263305 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1161 : END DO
1162 263305 : 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 263305 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1166 : END DO
1167 83874 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1168 366 : CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
1169 2342 : 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 2342 : alpha_scalar=1.0_dp, beta_scalar=qs_ot_env(1)%c_diis(i))
1173 : END DO
1174 2342 : 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 2342 : 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 40442 : IF (do_ener) THEN
1183 0 : DO ispin = 1, nspin
1184 0 : qs_ot_env(ispin)%ener_x(:) = 0.0_dp
1185 0 : DO i = 1, diis_bound
1186 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1187 0 : + qs_ot_env(1)%c_diis(i)*qs_ot_env(ispin)%ener_h_e(i, :)
1188 : END DO
1189 0 : DO i = 1, diis_bound
1190 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) &
1191 0 : + 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 40442 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1196 40442 : 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 40442 : tr_xold_gx = 0.0_dp
1201 40442 : tr_xnew_gx = 0.0_dp
1202 40442 : IF (do_ks) THEN
1203 83874 : DO ispin = 1, nspin
1204 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_h_x(j)%matrix, &
1205 43432 : qs_ot_env(ispin)%matrix_gx, tmp)
1206 43432 : tr_xold_gx = tr_xold_gx + tmp
1207 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, &
1208 43432 : qs_ot_env(ispin)%matrix_gx, tmp)
1209 43432 : tr_xnew_gx = tr_xnew_gx + tmp
1210 83874 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1211 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_h_x(j)%matrix, &
1212 366 : qs_ot_env(ispin)%rot_mat_gx, tmp)
1213 366 : tr_xold_gx = tr_xold_gx + 0.5_dp*tmp
1214 : CALL dbcsr_dot(qs_ot_env(ispin)%rot_mat_x, &
1215 366 : qs_ot_env(ispin)%rot_mat_gx, tmp)
1216 366 : tr_xnew_gx = tr_xnew_gx + 0.5_dp*tmp
1217 : END IF
1218 : END DO
1219 : END IF
1220 40442 : IF (do_ener) THEN
1221 0 : DO ispin = 1, nspin
1222 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_h_x(j, :), qs_ot_env(ispin)%ener_gx(:))
1223 0 : tr_xold_gx = tr_xold_gx + tmp
1224 0 : tmp = DOT_PRODUCT(qs_ot_env(ispin)%ener_x(:), qs_ot_env(ispin)%ener_gx(:))
1225 0 : tr_xnew_gx = tr_xnew_gx + tmp
1226 : END DO
1227 : END IF
1228 40442 : overlap = (tr_xnew_gx - tr_xold_gx)
1229 : ! OK, bad luck, take a SD step along the preconditioned gradient
1230 40442 : 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 840 : qs_ot_env(1)%OT_METHOD_FULL = "OT SD"
1238 840 : IF (do_ks) THEN
1239 1794 : DO ispin = 1, nspin
1240 954 : 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 954 : 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 954 : 1.0_dp, 1.0_dp)
1247 1794 : IF (qs_ot_env(ispin)%settings%do_rotation) THEN
1248 28 : 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 28 : 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 28 : 1.0_dp, 1.0_dp)
1255 : END IF
1256 : END DO
1257 : END IF
1258 840 : IF (do_ener) THEN
1259 0 : DO ispin = 1, nspin
1260 0 : qs_ot_env(ispin)%ener_x(:) = 0._dp
1261 0 : qs_ot_env(ispin)%ener_x(:) = qs_ot_env(ispin)%ener_x(:) + qs_ot_env(ispin)%ener_h_e(j, :)
1262 0 : 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 40444 : CALL timestop(handle)
1268 :
1269 40444 : 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 176 : 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 176 : 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 176 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f, x
1287 176 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: G, S
1288 : TYPE(cp_logger_type), POINTER :: logger
1289 :
1290 176 : eta = qs_ot_env(1)%settings%broyden_eta
1291 176 : omega = qs_ot_env(1)%settings%broyden_omega
1292 176 : sigma_dec = qs_ot_env(1)%settings%broyden_sigma_decrease
1293 176 : sigma_min = qs_ot_env(1)%settings%broyden_sigma_min
1294 176 : forget_history = qs_ot_env(1)%settings%broyden_forget_history
1295 176 : adaptive_sigma = qs_ot_env(1)%settings%broyden_adaptive_sigma
1296 176 : enable_flip = qs_ot_env(1)%settings%broyden_enable_flip
1297 :
1298 176 : do_ks = qs_ot_env(1)%settings%ks
1299 176 : do_ener = qs_ot_env(1)%settings%do_ener
1300 :
1301 176 : beta = qs_ot_env(1)%settings%broyden_beta
1302 176 : gamma = qs_ot_env(1)%settings%broyden_gamma
1303 176 : IF (adaptive_sigma) THEN
1304 176 : IF (qs_ot_env(1)%broyden_adaptive_sigma .LT. 0.0_dp) THEN
1305 6 : 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 0 : sigma = qs_ot_env(1)%settings%broyden_sigma
1311 : END IF
1312 :
1313 : ! simplify our life....
1314 176 : IF (do_ener .OR. .NOT. do_ks .OR. qs_ot_env(1)%settings%do_rotation) THEN
1315 0 : CPABORT("Not yet implemented")
1316 : END IF
1317 : !
1318 176 : nspin = SIZE(qs_ot_env)
1319 :
1320 176 : diis_m = qs_ot_env(1)%settings%diis_m
1321 :
1322 176 : IF (qs_ot_env(1)%diis_iter .LT. diis_m) THEN
1323 84 : diis_bound = qs_ot_env(1)%diis_iter + 1
1324 : ELSE
1325 92 : diis_bound = diis_m
1326 : END IF
1327 :
1328 : ! We want x:s, f:s and one random vector
1329 176 : k = 2*diis_bound + 1
1330 704 : ALLOCATE (S(k, k))
1331 528 : ALLOCATE (G(k, k))
1332 528 : ALLOCATE (f(k))
1333 352 : ALLOCATE (x(k))
1334 528 : ALLOCATE (circ_index(diis_bound))
1335 30400 : G = 0.0
1336 2272 : DO i = 1, k
1337 2272 : G(i, i) = sigma
1338 : END DO
1339 30400 : S = 0.0
1340 :
1341 176 : j = MOD(qs_ot_env(1)%diis_iter, diis_m) + 1 ! index in the circular array
1342 :
1343 352 : DO ispin = 1, nspin
1344 352 : CALL dbcsr_copy(qs_ot_env(ispin)%matrix_h_x(j)%matrix, qs_ot_env(ispin)%matrix_x)
1345 : END DO
1346 :
1347 176 : IF (ASSOCIATED(qs_ot_env(1)%preconditioner)) THEN
1348 176 : qs_ot_env(1)%gnorm = 0.0_dp
1349 352 : DO ispin = 1, nspin
1350 : CALL apply_preconditioner(qs_ot_env(ispin)%preconditioner, &
1351 176 : 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 176 : tmp)
1354 352 : qs_ot_env(1)%gnorm = qs_ot_env(1)%gnorm + tmp
1355 : END DO
1356 176 : IF (qs_ot_env(1)%gnorm .LT. 0.0_dp) THEN
1357 0 : WRITE (cp_logger_get_default_unit_nr(logger), *) "WARNING Preconditioner not positive definite !"
1358 : END IF
1359 352 : DO ispin = 1, nspin
1360 352 : CALL dbcsr_scale(qs_ot_env(ispin)%matrix_h_e(j)%matrix, -1.0_dp)
1361 : END DO
1362 : ELSE
1363 0 : qs_ot_env(1)%gnorm = 0.0_dp
1364 0 : DO ispin = 1, nspin
1365 0 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, qs_ot_env(ispin)%matrix_gx, tmp)
1366 0 : 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 0 : qs_ot_env(ispin)%matrix_gx, alpha_scalar=0.0_dp, beta_scalar=-1.0_dp)
1369 : END DO
1370 : END IF
1371 :
1372 176 : k = 0
1373 : n = 0
1374 176 : CALL dbcsr_get_info(qs_ot_env(1)%matrix_x, nfullrows_total=n)
1375 352 : DO ispin = 1, nspin
1376 176 : CALL dbcsr_get_info(qs_ot_env(ispin)%matrix_x, nfullcols_total=itmp)
1377 352 : k = k + itmp
1378 : END DO
1379 :
1380 : ! Handling the case of no free variables to optimize
1381 176 : IF (INT(n, KIND=int_8)*INT(k, KIND=int_8) /= 0) THEN
1382 176 : qs_ot_env(1)%delta = SQRT(ABS(qs_ot_env(1)%gnorm)/(INT(n, KIND=int_8)*INT(k, KIND=int_8)))
1383 176 : qs_ot_env(1)%gradient = -qs_ot_env(1)%gnorm
1384 : ELSE
1385 0 : qs_ot_env(1)%delta = 0.0_dp
1386 0 : qs_ot_env(1)%gradient = 0.0_dp
1387 : END IF
1388 :
1389 176 : IF (diis_bound == diis_m) THEN
1390 816 : DO i = 1, diis_bound
1391 816 : circ_index(i) = MOD(j + i - 1, diis_m) + 1
1392 : END DO
1393 : ELSE
1394 320 : DO i = 1, diis_bound
1395 320 : circ_index(i) = i
1396 : END DO
1397 : END IF
1398 :
1399 30400 : S = 0.0_dp
1400 352 : DO ispin = 1, nspin
1401 176 : CALL dbcsr_init_random(qs_ot_env(ispin)%matrix_x)
1402 1136 : DO i = 1, diis_bound
1403 : CALL dbcsr_dot( &
1404 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1405 960 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, tmp)
1406 960 : S(i, i) = S(i, i) + tmp
1407 : CALL dbcsr_dot( &
1408 : qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, &
1409 960 : qs_ot_env(ispin)%matrix_h_e(circ_index(i))%matrix, tmp)
1410 960 : 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 960 : qs_ot_env(ispin)%matrix_x, tmp)
1414 960 : S(i, 2*diis_bound + 1) = S(i, 2*diis_bound + 1) + tmp
1415 960 : 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 960 : qs_ot_env(ispin)%matrix_x, tmp)
1419 960 : S(i + diis_bound, 2*diis_bound + 1) = S(i + diis_bound, 2*diis_bound + 1) + tmp
1420 960 : S(i + diis_bound, 2*diis_bound + 1) = S(2*diis_bound + 1, diis_bound + i)
1421 3494 : 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 2534 : tmp)
1426 2534 : S(i, k) = S(i, k) + tmp
1427 2534 : 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 2534 : tmp)
1432 2534 : S(diis_bound + i, diis_bound + k) = S(diis_bound + i, diis_bound + k) + tmp
1433 3494 : S(diis_bound + k, diis_bound + i) = S(diis_bound + i, diis_bound + k)
1434 : END DO
1435 8124 : DO k = 1, diis_bound
1436 : CALL dbcsr_dot( &
1437 : qs_ot_env(ispin)%matrix_h_x(circ_index(i))%matrix, &
1438 6028 : qs_ot_env(ispin)%matrix_h_e(circ_index(k))%matrix, tmp)
1439 6028 : S(i, k + diis_bound) = S(i, k + diis_bound) + tmp
1440 6988 : S(k + diis_bound, i) = S(i, k + diis_bound)
1441 : END DO
1442 : END DO
1443 176 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_x, tmp)
1444 352 : 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 176 : k = 2*diis_bound + 1
1449 176 : tmp = SQRT(S(k, k))
1450 2272 : S(k, :) = S(k, :)/tmp
1451 2272 : S(:, k) = S(:, k)/tmp
1452 :
1453 176 : IF (diis_bound .GT. 1) THEN
1454 162 : tmp = 0.0_dp
1455 162 : tmp2 = 0.0_dp
1456 162 : i = diis_bound
1457 324 : 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 162 : tmp)
1463 162 : 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 162 : tmp)
1468 162 : 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 162 : tmp)
1473 162 : 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 162 : tmp)
1478 324 : tmp2 = tmp2 + tmp
1479 : END DO
1480 162 : qs_ot_env(1)%c_broy(i - 1) = tmp2
1481 : END IF
1482 :
1483 176 : qs_ot_env(1)%energy_h(j) = qs_ot_env(1)%etotal
1484 :
1485 : ! If we went uphill, do backtracking line search
1486 1136 : i = MINLOC(qs_ot_env(1)%energy_h(1:diis_bound), dim=1)
1487 176 : IF (i .NE. j) THEN
1488 26 : sigma = sigma_dec*sigma
1489 26 : qs_ot_env(1)%OT_METHOD_FULL = "OT BTRK"
1490 52 : DO ispin = 1, nspin
1491 26 : 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 26 : 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 52 : alpha_scalar=1.0_dp, beta_scalar=gamma)
1498 : END DO
1499 : ELSE
1500 : ! Construct G
1501 778 : DO i = 2, diis_bound
1502 9208 : f = 0.0
1503 9208 : x = 0.0
1504 : ! f is df_i
1505 628 : x(i) = 1.0
1506 628 : x(i - 1) = -1.0
1507 : ! x is dx_i
1508 628 : f(diis_bound + i) = 1.0
1509 628 : f(diis_bound + i - 1) = -1.0
1510 628 : tmp = 1.0_dp
1511 : ! We want a pos def Hessian
1512 628 : IF (enable_flip) THEN
1513 628 : IF (qs_ot_env(1)%c_broy(i - 1) .GT. 0) THEN
1514 : !qs_ot_env(1)%OT_METHOD_FULL="OT FLIP"
1515 2 : tmp = -1.0_dp
1516 : END IF
1517 : END IF
1518 :
1519 : ! get dx-Gdf
1520 391524 : x(:) = tmp*x - MATMUL(G, f)
1521 : ! dfSdf
1522 : ! we calculate matmul(S, f) twice. They're small...
1523 391524 : tmp = DOT_PRODUCT(f, MATMUL(S, f))
1524 : ! NOTE THAT S IS SYMMETRIC !!!
1525 391524 : f(:) = MATMUL(S, f)/tmp
1526 : ! the spread is an outer vector product
1527 130658 : G(:, :) = G + SPREAD(x, dim=2, ncopies=SIZE(f))*SPREAD(f, dim=1, ncopies=SIZE(x))
1528 : END DO
1529 1856 : f = 0.0_dp
1530 150 : f(2*diis_bound) = 1.0_dp
1531 94380 : x(:) = -beta*MATMUL(G, f)
1532 :
1533 : ! OK, add the vectors now, this sums up to the proposed step
1534 300 : DO ispin = 1, nspin
1535 150 : CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
1536 928 : 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 928 : alpha_scalar=1.0_dp, beta_scalar=-x(i + diis_bound))
1540 : END DO
1541 1078 : 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 928 : alpha_scalar=1.0_dp, beta_scalar=x(i))
1545 : END DO
1546 : END DO
1547 :
1548 150 : IF (adaptive_sigma) THEN
1549 150 : tmp = new_sigma(G, S, diis_bound)
1550 : !tmp = tmp * qs_ot_env ( 1 ) % settings % broyden_sigma
1551 150 : tmp = tmp*eta
1552 150 : sigma = MIN(omega*sigma, tmp)
1553 : END IF
1554 :
1555 : ! compute the inner product of direction of the step and gradient
1556 150 : tmp = 0.0_dp
1557 300 : DO ispin = 1, nspin
1558 : CALL dbcsr_dot(qs_ot_env(ispin)%matrix_gx, &
1559 : qs_ot_env(ispin)%matrix_x, &
1560 150 : tmp2)
1561 300 : tmp = tmp + tmp2
1562 : END DO
1563 :
1564 300 : DO ispin = 1, nspin
1565 : ! if the direction of the step is not in direction of the gradient,
1566 : ! change step sign
1567 300 : IF (tmp .GE. 0.0_dp) THEN
1568 8 : qs_ot_env(1)%OT_METHOD_FULL = "OT TURN"
1569 8 : IF (forget_history) THEN
1570 0 : qs_ot_env(1)%diis_iter = 0
1571 : END IF
1572 8 : 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 8 : 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 142 : 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 176 : DEALLOCATE (S, G, f, x, circ_index)
1586 :
1587 : ! update for next round
1588 176 : qs_ot_env(1)%diis_iter = qs_ot_env(1)%diis_iter + 1
1589 176 : qs_ot_env(1)%broyden_adaptive_sigma = MAX(sigma, sigma_min)
1590 :
1591 352 : END SUBROUTINE ot_broyden_step
1592 :
1593 : ! **************************************************************************************************
1594 : !> \brief ...
1595 : !> \param G ...
1596 : !> \param S ...
1597 : !> \param n ...
1598 : !> \return ...
1599 : ! **************************************************************************************************
1600 150 : 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 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H
1612 :
1613 600 : ALLOCATE (H(n, n))
1614 150 : CALL hess_G(G, S, H, n)
1615 450 : ALLOCATE (eigv(n))
1616 150 : 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 1836 : 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 150 : DEALLOCATE (H, eigv)
1631 150 : END FUNCTION new_sigma
1632 :
1633 : ! **************************************************************************************************
1634 : !> \brief ...
1635 : !> \param G ...
1636 : !> \param S ...
1637 : !> \param H ...
1638 : !> \param n ...
1639 : ! **************************************************************************************************
1640 150 : 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 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: v
1653 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Q
1654 :
1655 150 : i = SIZE(G, 1)
1656 150 : k = SIZE(H, 1)
1657 600 : ALLOCATE (Q(i, k))
1658 450 : ALLOCATE (v(i))
1659 5682 : H = 0.0_dp
1660 11214 : Q = 0.0_dp
1661 :
1662 1856 : Q(:, 1) = 1.0_dp
1663 27846 : tmp = SQRT(DOT_PRODUCT(Q(:, 1), MATMUL(S, Q(:, 1))))
1664 11214 : Q(:, :) = Q(:, :)/tmp
1665 :
1666 914 : DO i = 1, k
1667 162368 : v(:) = MATMUL(G, Q(:, i))
1668 3476 : DO j = 1, i
1669 658356 : H(j, i) = DOT_PRODUCT(Q(:, j), MATMUL(S, v))
1670 41072 : v(:) = v - H(j, i)*Q(:, j)
1671 : END DO
1672 914 : IF (i .LT. k) THEN
1673 147286 : tmp = DOT_PRODUCT(v, MATMUL(S, v))
1674 622 : IF (tmp .LE. 0.0_dp) THEN
1675 0 : n = i
1676 0 : EXIT
1677 : END IF
1678 622 : tmp = SQRT(tmp)
1679 : ! Lucky breakdown
1680 622 : IF (ABS(tmp) .LT. 1e-9_dp) THEN
1681 4 : n = i
1682 4 : EXIT
1683 : END IF
1684 618 : H(i + 1, i) = tmp
1685 9048 : Q(:, i + 1) = v/H(i + 1, i)
1686 : END IF
1687 : END DO
1688 :
1689 150 : DEALLOCATE (Q, v)
1690 150 : END SUBROUTINE hess_G
1691 :
1692 5364 : END MODULE qs_ot_minimizer
|
