(git:34ef472)
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 ! **************************************************************************************************
12 MODULE 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 
114 CONTAINS
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 
723 END MODULE linesearch
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