(git:374b731)
Loading...
Searching...
No Matches
linesearch.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 A generic framework to calculate step lengths for 1D line search
10!> \author Ole Schuett
11! **************************************************************************************************
12MODULE linesearch
13 USE cp_log_handling, ONLY: cp_get_default_logger,&
14 cp_logger_type
15 USE cp_output_handling, ONLY: add_last_numeric,&
16 cp_print_key_section_create,&
17 cp_print_key_unit_nr,&
18 low_print_level
19 USE input_keyword_types, ONLY: keyword_create,&
20 keyword_release,&
21 keyword_type
22 USE input_section_types, ONLY: section_add_keyword,&
23 section_add_subsection,&
24 section_create,&
25 section_release,&
26 section_type,&
27 section_vals_type,&
28 section_vals_val_get
29 USE kinds, ONLY: dp
30 USE string_utilities, ONLY: s2a
31#include "./base/base_uses.f90"
32
33 IMPLICIT NONE
34
35 PRIVATE
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'linesearch'
38
39 PUBLIC :: linesearch_type
40
41 INTEGER, PARAMETER :: linesearch_method_adapt = 1
42 INTEGER, PARAMETER :: linesearch_method_2pnt = 2
43 INTEGER, PARAMETER :: linesearch_method_3pnt = 3
44 INTEGER, PARAMETER :: linesearch_method_gold = 4
45 INTEGER, PARAMETER :: linesearch_method_none = 5
46
47 TYPE linesearch_2pnt_type
48 REAL(KIND=dp), DIMENSION(2) :: energies = 0.0
49 REAL(KIND=dp) :: scan_step = 0.0
50 REAL(KIND=dp) :: last_step_size = 0.0
51 REAL(KIND=dp) :: max_step_size = 0.0
52 INTEGER :: count = 1
53 END TYPE linesearch_2pnt_type
54
55 TYPE linesearch_3pnt_type
56 REAL(KIND=dp), DIMENSION(3) :: energies = 0.0
57 REAL(KIND=dp), DIMENSION(3) :: scan_steps = 0.0
58 REAL(KIND=dp) :: last_step_size = 0.0
59 REAL(KIND=dp) :: max_step_size = 0.0
60 REAL(KIND=dp) :: tiny_step_size = 0.0
61 INTEGER :: count = 1
62 END TYPE linesearch_3pnt_type
63
64 TYPE linesearch_adapt_type
65 REAL(KIND=dp) :: last_step_size = 0.0
66 REAL(KIND=dp) :: left_x = 0.0
67 REAL(KIND=dp) :: middle_x = 0.0
68 REAL(KIND=dp) :: right_x = 0.0
69 REAL(KIND=dp) :: left_e = 0.0
70 REAL(KIND=dp) :: middle_e = 0.0
71 REAL(KIND=dp) :: right_e = 0.0
72 LOGICAL :: have_left = .false.
73 LOGICAL :: have_middle = .false.
74 LOGICAL :: have_right = .false.
75 INTEGER :: count = 0
76 END TYPE linesearch_adapt_type
77
78 TYPE linesearch_gold_type
79 REAL(KIND=dp) :: scan_steps = 0.0
80 REAL(KIND=dp) :: last_step_size = 0.0
81 REAL(KIND=dp) :: eps_step_size = 0.0
82 REAL(KIND=dp) :: left_x = 0.0
83 REAL(KIND=dp) :: middle_x = 0.0
84 REAL(KIND=dp) :: right_x = 0.0
85 REAL(KIND=dp) :: left_e = 0.0
86 REAL(KIND=dp) :: middle_e = 0.0
87 REAL(KIND=dp) :: right_e = 0.0
88 LOGICAL :: have_left = .false.
89 LOGICAL :: have_middle = .false.
90 LOGICAL :: have_right = .false.
91 LOGICAL :: gave_up = .false.
92 END TYPE linesearch_gold_type
93
94 TYPE linesearch_type
95 PRIVATE
96 REAL(kind=dp), PUBLIC :: step_size = 0.0_dp
97 LOGICAL, PUBLIC :: starts = .false.
98 TYPE(linesearch_adapt_type), POINTER :: state_adapt => null()
99 TYPE(linesearch_2pnt_type), POINTER :: state_2pnt => null()
100 TYPE(linesearch_3pnt_type), POINTER :: state_3pnt => null()
101 TYPE(linesearch_gold_type), POINTER :: state_gold => null()
102 INTEGER :: iw = -1
103 INTEGER :: method = -1
104 CHARACTER(LEN=10) :: label = ""
105 REAL(kind=dp) :: init_step_size = 0.0_dp
106 REAL(dp) :: eps_step_size = 0.0_dp
107 REAL(dp) :: max_step_size = 0.0_dp
108 REAL(dp) :: tiny_step_size = 0.0_dp
109 END TYPE linesearch_type
110
111 PUBLIC :: linesearch_create_section, linesearch_init, linesearch_finalize
112 PUBLIC :: linesearch_step, linesearch_reset
113
114CONTAINS
115
116! **************************************************************************************************
117!> \brief Declare the line search input section.
118!> \param section ...
119!> \author Ole Schuett
120! **************************************************************************************************
121 SUBROUTINE linesearch_create_section(section)
122 TYPE(section_type), POINTER :: section
123
124 TYPE(keyword_type), POINTER :: keyword
125 TYPE(section_type), POINTER :: print_section, printkey
126
127 NULLIFY (keyword, print_section, printkey)
128
129 cpassert(.NOT. ASSOCIATED(section))
130 CALL section_create(section, __location__, name="LINE_SEARCH", repeats=.false., &
131 description="Detail settings or linesearch method.")
132
133 CALL keyword_create( &
134 keyword, __location__, name="METHOD", &
135 description="Linesearch method.", &
136 default_i_val=linesearch_method_adapt, &
137 enum_c_vals=s2a("ADAPT", "3PNT", "2PNT", "GOLD", "NONE"), &
138 enum_desc=s2a("extrapolates usually based on 3 points, uses additional points on demand, very robust.", &
139 "extrapolate based on 3 points", &
140 "extrapolate based on 2 points and the slope (super fast, but might get stuck at saddle points)", &
141 "perform 1D golden section search of the minimum (very expensive)", &
142 "always take steps of fixed INITIAL_STEP_SIZE"), &
143 enum_i_vals=(/linesearch_method_adapt, linesearch_method_3pnt, linesearch_method_2pnt, &
144 linesearch_method_gold, linesearch_method_none/))
145 CALL section_add_keyword(section, keyword)
146 CALL keyword_release(keyword)
147
148 CALL keyword_create(keyword, __location__, name="INITIAL_STEP_SIZE", &
149 description="Initial step length", &
150 default_r_val=0.1_dp)
151 CALL section_add_keyword(section, keyword)
152 CALL keyword_release(keyword)
153
154 CALL keyword_create(keyword, __location__, name="MAX_STEP_SIZE", &
155 description="Maximum step length", &
156 default_r_val=3.0_dp)
157 CALL section_add_keyword(section, keyword)
158 CALL keyword_release(keyword)
159
160 CALL keyword_create(keyword, __location__, name="TINY_STEP_SIZE", &
161 description="Step length taken if negative step is suggested.", &
162 default_r_val=0.002_dp)
163 CALL section_add_keyword(section, keyword)
164 CALL keyword_release(keyword)
165
166 CALL keyword_create(keyword, __location__, name="EPS_STEP_SIZE", &
167 description="Convergence criterion of GOLD method.", &
168 default_r_val=0.1_dp)
169 CALL section_add_keyword(section, keyword)
170 CALL keyword_release(keyword)
171
172 CALL section_create(print_section, __location__, name="PRINT", &
173 description="Print section", &
174 n_keywords=0, n_subsections=1, repeats=.true.)
175
176 CALL cp_print_key_section_create(printkey, __location__, "RUN_INFO", &
177 description="General run information", &
178 print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
179
180 CALL section_add_subsection(print_section, printkey)
181 CALL section_release(printkey)
182 CALL section_add_subsection(section, print_section)
183 CALL section_release(print_section)
184
185 END SUBROUTINE linesearch_create_section
186
187! **************************************************************************************************
188!> \brief Initialize linesearch from given input section
189!> \param this ...
190!> \param section ...
191!> \param label ...
192!> \author Ole Schuett
193! **************************************************************************************************
194 SUBROUTINE linesearch_init(this, section, label)
195 TYPE(linesearch_type), INTENT(INOUT) :: this
196 TYPE(section_vals_type), POINTER :: section
197 CHARACTER(LEN=*) :: label
198
199 TYPE(cp_logger_type), POINTER :: logger
200
201 CALL section_vals_val_get(section, "METHOD", i_val=this%method)
202 CALL section_vals_val_get(section, "INITIAL_STEP_SIZE", r_val=this%init_step_size)
203 CALL section_vals_val_get(section, "MAX_STEP_SIZE", r_val=this%max_step_size)
204 CALL section_vals_val_get(section, "TINY_STEP_SIZE", r_val=this%tiny_step_size)
205 CALL section_vals_val_get(section, "EPS_STEP_SIZE", r_val=this%eps_step_size)
206
207 cpassert(len_trim(label) <= 10)
208 this%label = label
209 logger => cp_get_default_logger()
210 this%iw = cp_print_key_unit_nr(logger, section, "PRINT%RUN_INFO", &
211 extension=".linesearchlog")
212
213 CALL linesearch_init_low(this)
214
215 END SUBROUTINE linesearch_init
216
217! **************************************************************************************************
218!> \brief Helper routine to (re)-initialize line search machinery
219!> \param this ...
220!> \author Ole Schuett
221! **************************************************************************************************
222 SUBROUTINE linesearch_init_low(this)
223 TYPE(linesearch_type), INTENT(INOUT) :: this
224
225 this%step_size = 0.0_dp
226 this%starts = .true.
227
228 SELECT CASE (this%method)
229 CASE (linesearch_method_adapt)
230 ALLOCATE (this%state_adapt)
231 this%state_adapt%last_step_size = this%init_step_size
232 CASE (linesearch_method_2pnt)
233 ALLOCATE (this%state_2pnt)
234 this%state_2pnt%max_step_size = this%max_step_size
235 this%state_2pnt%last_step_size = this%init_step_size
236 CASE (linesearch_method_3pnt)
237 ALLOCATE (this%state_3pnt)
238 this%state_3pnt%last_step_size = this%init_step_size
239 this%state_3pnt%max_step_size = this%max_step_size
240 this%state_3pnt%tiny_step_size = this%tiny_step_size
241 CASE (linesearch_method_gold)
242 ALLOCATE (this%state_gold)
243 this%state_gold%last_step_size = this%init_step_size
244 this%state_gold%eps_step_size = this%eps_step_size
245 CASE (linesearch_method_none)
246 ! nothing todo
247 CASE DEFAULT
248 cpabort("unknown method")
249 END SELECT
250
251 END SUBROUTINE linesearch_init_low
252
253! **************************************************************************************************
254!> \brief Finzalize line search machinery
255!> \param this ...
256!> \author Ole Schuett
257! **************************************************************************************************
258 SUBROUTINE linesearch_finalize(this)
259 TYPE(linesearch_type), INTENT(INOUT) :: this
260
261 SELECT CASE (this%method)
262 CASE (linesearch_method_adapt)
263 DEALLOCATE (this%state_adapt)
264 CASE (linesearch_method_2pnt)
265 DEALLOCATE (this%state_2pnt)
266 CASE (linesearch_method_3pnt)
267 DEALLOCATE (this%state_3pnt)
268 CASE (linesearch_method_gold)
269 DEALLOCATE (this%state_gold)
270 CASE (linesearch_method_none)
271 ! nothing todo
272 CASE DEFAULT
273 cpabort("unknown method")
274 END SELECT
275
276 !TODO: should finish printkey, but don't have the section here
277 END SUBROUTINE linesearch_finalize
278
279! **************************************************************************************************
280!> \brief Reset line search to initial state
281!> \param this ...
282!> \author Ole Schuett
283! **************************************************************************************************
284 SUBROUTINE linesearch_reset(this)
285 TYPE(linesearch_type), INTENT(INOUT) :: this
286
287 CALL linesearch_finalize(this)
288 CALL linesearch_init_low(this)
289 END SUBROUTINE linesearch_reset
290
291! **************************************************************************************************
292!> \brief Calculate step length of next line search step.
293!> \param this ...
294!> \param energy ...
295!> \param slope ...
296!> \author Ole Schuett
297! **************************************************************************************************
298 SUBROUTINE linesearch_step(this, energy, slope)
299 TYPE(linesearch_type), INTENT(INOUT) :: this
300 REAL(kind=dp), INTENT(IN) :: energy, slope
301
302 LOGICAL :: is_done
303 REAL(kind=dp) :: step_size
304
305 SELECT CASE (this%method)
306 CASE (linesearch_method_adapt)
307 CALL linesearch_adapt(this%state_adapt, energy, step_size, is_done, this%iw, trim(this%label))
308 CASE (linesearch_method_2pnt)
309 CALL linesearch_2pnt(this%state_2pnt, energy, slope, step_size, is_done, this%iw, trim(this%label))
310 CASE (linesearch_method_3pnt)
311 CALL linesearch_3pnt(this%state_3pnt, energy, step_size, is_done, this%iw, trim(this%label))
312 CASE (linesearch_method_gold)
313 CALL linesearch_gold(this%state_gold, energy, step_size, is_done, this%iw, trim(this%label))
314 CASE (linesearch_method_none)
315 step_size = this%init_step_size ! take steps of fixed length
316 is_done = .true.
317 CASE DEFAULT
318 cpabort("unknown method")
319 END SELECT
320
321 this%step_size = step_size
322 this%starts = is_done
323 END SUBROUTINE linesearch_step
324
325! **************************************************************************************************
326!> \brief Perform a 2pnt linesearch
327!> \param this ...
328!> \param energy ...
329!> \param slope ...
330!> \param step_size ...
331!> \param is_done ...
332!> \param unit_nr ...
333!> \param label ...
334!> \author Ole Schuett
335! **************************************************************************************************
336 SUBROUTINE linesearch_2pnt(this, energy, slope, step_size, is_done, unit_nr, label)
337 TYPE(linesearch_2pnt_type), INTENT(INOUT) :: this
338 REAL(kind=dp), INTENT(IN) :: energy, slope
339 REAL(kind=dp), INTENT(OUT) :: step_size
340 LOGICAL, INTENT(OUT) :: is_done
341 INTEGER, INTENT(IN) :: unit_nr
342 CHARACTER(len=*), INTENT(IN) :: label
343
344 REAL(kind=dp) :: a, b, c, pred_energy, x2
345
346 this%energies(this%count) = energy
347 is_done = .false.
348
349 SELECT CASE (this%count)
350 CASE (1)
351 step_size = 2.0_dp*this%last_step_size
352 this%scan_step = step_size
353 this%count = 2
354 CASE (2)
355 c = this%energies(1)
356 b = -slope
357 x2 = this%scan_step
358 a = (this%energies(2) - b*x2 - c)/(x2**2)
359
360 IF (a < 0.0_dp) THEN
361 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| had to quench curvature"
362 a = 1.0e-15_dp
363 END IF
364
365 step_size = -b/(2.0_dp*a)
366 pred_energy = a*step_size**2 + b*step_size + c
367
368 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 2pnt suggested step_size: ", step_size
369 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 2pnt predicted energy", pred_energy
370
371 IF (pred_energy > this%energies(1) .OR. pred_energy > this%energies(2)) THEN
372 cpabort(label//"LS| predicted energy not below test points")
373 END IF
374
375 IF (step_size > this%max_step_size) THEN
376 step_size = this%max_step_size
377 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| limiting step_size to MAX_STEP_SIZE"
378 END IF
379
380 this%last_step_size = step_size
381 this%count = 1
382 is_done = .true.
383 CASE DEFAULT
384 cpabort("this should not happen")
385 END SELECT
386
387 END SUBROUTINE linesearch_2pnt
388
389! **************************************************************************************************
390!> \brief Perform a 3pnt linesearch
391!> \param this ...
392!> \param energy ...
393!> \param step_size ...
394!> \param is_done ...
395!> \param unit_nr ...
396!> \param label ...
397!> \author Ole Schuett
398! **************************************************************************************************
399 SUBROUTINE linesearch_3pnt(this, energy, step_size, is_done, unit_nr, label)
400 TYPE(linesearch_3pnt_type), INTENT(INOUT) :: this
401 REAL(kind=dp), INTENT(IN) :: energy
402 REAL(kind=dp), INTENT(OUT) :: step_size
403 LOGICAL, INTENT(OUT) :: is_done
404 INTEGER, INTENT(IN) :: unit_nr
405 CHARACTER(len=*), INTENT(IN) :: label
406
407 REAL(kind=dp) :: a, b, c, denom, pred_energy, x1, x2, x3, &
408 y1, y2, y3
409
410 this%energies(this%count) = energy
411 is_done = .false.
412
413 SELECT CASE (this%count)
414 CASE (1)
415 step_size = (2.0_dp/3.0_dp)*this%last_step_size
416 IF (step_size < this%tiny_step_size) THEN
417 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| initial step size too small, using TINY_STEP_SIZE"
418 step_size = this%tiny_step_size
419 END IF
420 this%scan_steps(1) = 0.0_dp
421 this%scan_steps(2) = step_size
422 this%count = 2
423 CASE (2)
424 IF (this%energies(1) > this%energies(2)) THEN
425 step_size = 2.0_dp*this%scan_steps(2)
426 ELSE
427 step_size = 0.5_dp*this%scan_steps(2)
428 END IF
429 this%scan_steps(3) = step_size
430 this%count = 3
431 CASE (3)
432 ! fitting y = a*x^2 + b*x + c
433 y1 = this%energies(1); y2 = this%energies(2); y3 = this%energies(3)
434 x1 = this%scan_steps(1); x2 = this%scan_steps(2); x3 = this%scan_steps(3)
435
436 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt scan_steps: ", this%scan_steps
437 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt energies: ", this%energies
438
439 ! Cramer's Rule
440 denom = (x1 - x2)*(x1 - x3)*(x2 - x3)
441 a = (x3*(y2 - y1) + x2*(y1 - y3) + x1*(y3 - y2))/denom
442 b = (x3**2*(y1 - y2) + x2**2*(y3 - y1) + x1**2*(y2 - y3))/denom
443 c = (x2*x3*(x2 - x3)*y1 + x3*x1*(x3 - x1)*y2 + x1*x2*(x1 - x2)*y3)/denom
444
445 step_size = -b/(2.0_dp*a)
446 pred_energy = a*step_size**2 + b*step_size + c
447 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt suggested step_size: ", step_size
448 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt predicted energy", pred_energy
449
450 IF (a < 0) THEN
451 step_size = -2.0_dp*step_size
452 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| inverting step size"
453 END IF
454
455 IF (step_size < 0) THEN
456 step_size = this%tiny_step_size
457 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| makeing a step of size TINY_STEP_SIZE"
458 END IF
459
460 IF (step_size > this%max_step_size) THEN
461 step_size = this%max_step_size
462 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| limiting step_size to MAX_STEP_SIZE"
463 END IF
464
465 this%last_step_size = step_size
466 this%count = 1
467 is_done = .true.
468 CASE DEFAULT
469 cpabort("this should not happen")
470 END SELECT
471
472 END SUBROUTINE linesearch_3pnt
473
474! **************************************************************************************************
475!> \brief Perform an adaptive linesearch
476!> \param this ...
477!> \param energy ...
478!> \param step_size ...
479!> \param is_done ...
480!> \param unit_nr ...
481!> \param label ...
482!> \author Ole Schuett
483! **************************************************************************************************
484 SUBROUTINE linesearch_adapt(this, energy, step_size, is_done, unit_nr, label)
485 TYPE(linesearch_adapt_type), INTENT(INOUT) :: this
486 REAL(kind=dp), INTENT(IN) :: energy
487 REAL(kind=dp), INTENT(OUT) :: step_size
488 LOGICAL, INTENT(OUT) :: is_done
489 INTEGER, INTENT(IN) :: unit_nr
490 CHARACTER(len=*), INTENT(IN) :: label
491
492 REAL(kind=dp), PARAMETER :: grow_factor = 2.0_dp, &
493 shrink_factor = 0.5_dp
494
495 REAL(kind=dp) :: a, b, c, denom, pred_energy, x1, x2, x3, &
496 y1, y2, y3
497
498 is_done = .false.
499 this%count = this%count + 1
500
501 IF (.NOT. this%have_left) THEN
502 this%left_x = 0.0_dp
503 this%left_e = energy
504 this%have_left = .true.
505 step_size = this%last_step_size
506
507 ELSE IF (.NOT. (this%have_middle .OR. this%have_right)) THEN
508 IF (energy < this%left_e) THEN
509 this%middle_e = energy
510 this%middle_x = this%last_step_size
511 this%have_middle = .true.
512 step_size = this%middle_x*grow_factor
513 ELSE
514 this%right_e = energy
515 this%right_x = this%last_step_size
516 this%have_right = .true.
517 step_size = this%right_x*shrink_factor
518 END IF
519
520 ELSE IF (.NOT. this%have_right) THEN
521 IF (energy < this%middle_e) THEN
522 this%middle_e = energy
523 this%middle_x = this%last_step_size
524 step_size = this%middle_x*grow_factor
525 ELSE
526 this%right_e = energy
527 this%right_x = this%last_step_size
528 this%have_right = .true.
529 END IF
530
531 ELSE IF (.NOT. this%have_middle) THEN
532 IF (energy > this%left_e) THEN
533 this%right_e = energy
534 this%right_x = this%last_step_size
535 step_size = this%right_x*shrink_factor
536 ELSE
537 this%middle_e = energy
538 this%middle_x = this%last_step_size
539 this%have_middle = .true.
540 END IF
541 END IF
542
543 IF (this%count > 3) THEN
544 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| Need extra step"
545 END IF
546 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: ", this%have_left, this%have_middle, this%have_right
547 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: scan_steps: ", this%left_x, this%middle_x, this%right_x
548 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: energies: ", this%left_e, this%middle_e, this%right_e
549
550 IF (this%have_left .AND. this%have_middle .AND. this%have_right) THEN
551 ! fitting y = a*x^2 + b*x + c
552 y1 = this%left_e; y2 = this%middle_e; y3 = this%right_e
553 x1 = this%left_x; x2 = this%middle_x; x3 = this%right_x
554
555 ! Cramer's rule
556 denom = (x1 - x2)*(x1 - x3)*(x2 - x3)
557 a = (x3*(y2 - y1) + x2*(y1 - y3) + x1*(y3 - y2))/denom
558 b = (x3**2*(y1 - y2) + x2**2*(y3 - y1) + x1**2*(y2 - y3))/denom
559 c = (x2*x3*(x2 - x3)*y1 + x3*x1*(x3 - x1)*y2 + x1*x2*(x1 - x2)*y3)/denom
560
561 IF (abs(a) /= 0.0_dp) THEN
562 step_size = -b/(2.0_dp*a)
563 ELSE
564 step_size = 0.0_dp
565 END IF
566
567 cpassert(step_size >= 0.0_dp)
568 pred_energy = a*step_size**2 + b*step_size + c
569 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: suggested step_size: ", step_size
570 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: predicted energy", pred_energy
571
572 ! reset
573 is_done = .true.
574 this%count = 0
575 this%have_left = .false.
576 this%have_middle = .false.
577 this%have_right = .false.
578 this%left_e = 0.0
579 this%middle_e = 0.0
580 this%right_e = 0.0
581 this%left_x = 0.0
582 this%middle_x = 0.0
583 this%right_x = 0.0
584 END IF
585
586 this%last_step_size = step_size
587 END SUBROUTINE linesearch_adapt
588
589! **************************************************************************************************
590!> \brief Perform a gold linesearch
591!> \param this ...
592!> \param energy ...
593!> \param step_size ...
594!> \param is_done ...
595!> \param unit_nr ...
596!> \param label ...
597!> \author Ole Schuett
598! **************************************************************************************************
599 SUBROUTINE linesearch_gold(this, energy, step_size, is_done, unit_nr, label)
600 TYPE(linesearch_gold_type), INTENT(INOUT) :: this
601 REAL(kind=dp), INTENT(IN) :: energy
602 REAL(kind=dp), INTENT(OUT) :: step_size
603 LOGICAL, INTENT(OUT) :: is_done
604 INTEGER, INTENT(IN) :: unit_nr
605 CHARACTER(len=*), INTENT(IN) :: label
606
607 REAL(kind=dp), PARAMETER :: phi = (1.0_dp + sqrt(5.0_dp))/2.0_dp
608
609 REAL(kind=dp) :: a, b, d
610
611 is_done = .false.
612
613 IF (this%gave_up) &
614 cpabort("had to give up, should not be called again")
615
616 IF (.NOT. this%have_left) THEN
617 this%left_x = 0.0_dp
618 this%left_e = energy
619 this%have_left = .true.
620 step_size = this%last_step_size
621
622 ELSE IF (.NOT. (this%have_middle .OR. this%have_right)) THEN
623 IF (energy < this%left_e) THEN
624 this%middle_e = energy
625 this%middle_x = this%scan_steps
626 this%have_middle = .true.
627 step_size = this%middle_x*phi
628 ELSE
629 this%right_e = energy
630 this%right_x = this%scan_steps
631 this%have_right = .true.
632 step_size = this%right_x/phi
633 END IF
634
635 ELSE IF (.NOT. this%have_right) THEN
636 IF (energy < this%middle_e) THEN
637 this%middle_e = energy
638 this%middle_x = this%scan_steps
639 step_size = this%middle_x*phi
640 ELSE
641 this%right_e = energy
642 this%right_x = this%scan_steps
643 this%have_right = .true.
644 END IF
645
646 ELSE IF (.NOT. this%have_middle) THEN
647 IF (energy > this%left_e) THEN
648 this%right_e = energy
649 this%right_x = this%scan_steps
650 step_size = this%right_x/phi
651 ELSE
652 this%middle_e = energy
653 this%middle_x = this%scan_steps
654 this%have_middle = .true.
655 END IF
656
657 ELSE !up and running
658 a = this%middle_x - this%left_x
659 b = this%right_x - this%middle_x
660 IF (energy < this%middle_e) THEN
661 IF (a < b) THEN
662 this%left_e = this%middle_e
663 this%left_x = this%middle_x
664 ELSE
665 this%right_e = this%middle_e
666 this%right_x = this%middle_x
667 END IF
668 this%middle_e = energy
669 this%middle_x = this%scan_steps
670 ELSE
671 IF (a < b) THEN
672 this%right_e = energy
673 this%right_x = this%scan_steps
674 ELSE
675 this%left_e = energy
676 this%left_x = this%scan_steps
677 END IF
678 END IF
679 END IF
680
681 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%have_left, this%have_middle, this%have_right
682 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%left_x, this%middle_x, this%right_x
683 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%left_e, this%middle_e, this%right_e
684
685 IF (this%have_left .AND. this%have_middle .AND. this%have_right) THEN
686 a = this%middle_x - this%left_x
687 b = this%right_x - this%middle_x
688 IF (abs(min(a, b)*phi - max(a, b)) > 1.0e-10) &
689 cpabort("golden-ratio gone")
690
691 IF (a < b) THEN
692 step_size = this%middle_x + a/phi
693 ELSE
694 step_size = this%middle_x - b/phi
695 END IF
696
697 d = abs(this%right_x - this%left_x)/(abs(this%middle_x) + abs(step_size))
698 IF (d < this%eps_step_size) THEN
699 step_size = this%middle_x
700 this%last_step_size = step_size
701 is_done = .true.
702
703 IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold done, step-size: ", step_size
704
705 this%have_left = .false.
706 this%have_middle = .false.
707 this%have_right = .false.
708 this%left_e = 0.0
709 this%middle_e = 0.0
710 this%right_e = 0.0
711 this%left_x = 0.0
712 this%middle_x = 0.0
713 this%right_x = 0.0
714 END IF
715
716 END IF
717
718 IF (step_size < 1e-10) cpabort("linesearch failed / done")
719
720 this%scan_steps = step_size
721 END SUBROUTINE linesearch_gold
722
723END MODULE linesearch
string_utilities::s2a
Definition string_utilities.F:48
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_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::add_last_numeric
integer, parameter, public add_last_numeric
Definition cp_output_handling.F:70
cp_output_handling::cp_print_key_section_create
subroutine, public cp_print_key_section_create(print_key_section, location, name, description, print_level, each_iter_names, each_iter_values, add_last, filename, common_iter_levels, citations, unit_str)
creates a print_key section
Definition cp_output_handling.F:174
input_keyword_types
represents keywords in an input
Definition input_keyword_types.F:14
input_keyword_types::keyword_release
subroutine, public keyword_release(keyword)
releases the given keyword (see doc/ReferenceCounting.html)
Definition input_keyword_types.F:426
input_keyword_types::keyword_create
subroutine, public keyword_create(keyword, location, name, description, usage, type_of_var, n_var, repeats, variants, default_val, default_l_val, default_r_val, default_lc_val, default_c_val, default_i_val, default_l_vals, default_r_vals, default_c_vals, default_i_vals, lone_keyword_val, lone_keyword_l_val, lone_keyword_r_val, lone_keyword_c_val, lone_keyword_i_val, lone_keyword_l_vals, lone_keyword_r_vals, lone_keyword_c_vals, lone_keyword_i_vals, enum_c_vals, enum_i_vals, enum, enum_strict, enum_desc, unit_str, citations, deprecation_notice, removed)
creates a keyword object
Definition input_keyword_types.F:155
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_create
subroutine, public section_create(section, location, name, description, n_keywords, n_subsections, repeats, citations)
creates a list of keywords
Definition input_section_types.F:156
input_section_types::section_add_keyword
subroutine, public section_add_keyword(section, keyword)
adds a keyword to the given section
Definition input_section_types.F:504
input_section_types::section_add_subsection
subroutine, public section_add_subsection(section, subsection)
adds a subsection to the given section
Definition input_section_types.F:560
input_section_types::section_release
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
Definition input_section_types.F:231
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
linesearch
A generic framework to calculate step lengths for 1D line search.
Definition linesearch.F:12
linesearch::linesearch_create_section
subroutine, public linesearch_create_section(section)
Declare the line search input section.
Definition linesearch.F:122
linesearch::linesearch_finalize
subroutine, public linesearch_finalize(this)
Finzalize line search machinery.
Definition linesearch.F:259
linesearch::linesearch_init
subroutine, public linesearch_init(this, section, label)
Initialize linesearch from given input section.
Definition linesearch.F:195
linesearch::linesearch_reset
subroutine, public linesearch_reset(this)
Reset line search to initial state.
Definition linesearch.F:285
linesearch::linesearch_step
subroutine, public linesearch_step(this, energy, slope)
Calculate step length of next line search step.
Definition linesearch.F:299
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
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
input_keyword_types::keyword_type
represent a keyword in the input
Definition input_keyword_types.F:88
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
linesearch::linesearch_type
Definition linesearch.F:94