Line data Source code
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 9190 : 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 9190 : NULLIFY (keyword, print_section, printkey)
128 :
129 9190 : CPASSERT(.NOT. ASSOCIATED(section))
130 : CALL section_create(section, __LOCATION__, name="LINE_SEARCH", repeats=.FALSE., &
131 9190 : 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 9190 : linesearch_method_gold, linesearch_method_none/))
145 9190 : CALL section_add_keyword(section, keyword)
146 9190 : CALL keyword_release(keyword)
147 :
148 : CALL keyword_create(keyword, __LOCATION__, name="INITIAL_STEP_SIZE", &
149 : description="Initial step length", &
150 9190 : default_r_val=0.1_dp)
151 9190 : CALL section_add_keyword(section, keyword)
152 9190 : CALL keyword_release(keyword)
153 :
154 : CALL keyword_create(keyword, __LOCATION__, name="MAX_STEP_SIZE", &
155 : description="Maximum step length", &
156 9190 : default_r_val=3.0_dp)
157 9190 : CALL section_add_keyword(section, keyword)
158 9190 : 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 9190 : default_r_val=0.002_dp)
163 9190 : CALL section_add_keyword(section, keyword)
164 9190 : CALL keyword_release(keyword)
165 :
166 : CALL keyword_create(keyword, __LOCATION__, name="EPS_STEP_SIZE", &
167 : description="Convergence criterion of GOLD method.", &
168 9190 : default_r_val=0.1_dp)
169 9190 : CALL section_add_keyword(section, keyword)
170 9190 : CALL keyword_release(keyword)
171 :
172 : CALL section_create(print_section, __LOCATION__, name="PRINT", &
173 : description="Print section", &
174 9190 : 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 9190 : print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
179 :
180 9190 : CALL section_add_subsection(print_section, printkey)
181 9190 : CALL section_release(printkey)
182 9190 : CALL section_add_subsection(section, print_section)
183 9190 : CALL section_release(print_section)
184 :
185 9190 : 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 280 : 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 280 : CALL section_vals_val_get(section, "METHOD", i_val=this%method)
202 280 : CALL section_vals_val_get(section, "INITIAL_STEP_SIZE", r_val=this%init_step_size)
203 280 : CALL section_vals_val_get(section, "MAX_STEP_SIZE", r_val=this%max_step_size)
204 280 : CALL section_vals_val_get(section, "TINY_STEP_SIZE", r_val=this%tiny_step_size)
205 280 : CALL section_vals_val_get(section, "EPS_STEP_SIZE", r_val=this%eps_step_size)
206 :
207 280 : CPASSERT(LEN_TRIM(label) <= 10)
208 280 : this%label = label
209 280 : logger => cp_get_default_logger()
210 : this%iw = cp_print_key_unit_nr(logger, section, "PRINT%RUN_INFO", &
211 280 : extension=".linesearchlog")
212 :
213 280 : CALL linesearch_init_low(this)
214 :
215 280 : 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 526 : SUBROUTINE linesearch_init_low(this)
223 : TYPE(linesearch_type), INTENT(INOUT) :: this
224 :
225 526 : this%step_size = 0.0_dp
226 526 : this%starts = .TRUE.
227 :
228 928 : SELECT CASE (this%method)
229 : CASE (linesearch_method_adapt)
230 402 : ALLOCATE (this%state_adapt)
231 402 : this%state_adapt%last_step_size = this%init_step_size
232 : CASE (linesearch_method_2pnt)
233 64 : ALLOCATE (this%state_2pnt)
234 16 : this%state_2pnt%max_step_size = this%max_step_size
235 16 : this%state_2pnt%last_step_size = this%init_step_size
236 : CASE (linesearch_method_3pnt)
237 540 : ALLOCATE (this%state_3pnt)
238 60 : this%state_3pnt%last_step_size = this%init_step_size
239 60 : this%state_3pnt%max_step_size = this%max_step_size
240 60 : this%state_3pnt%tiny_step_size = this%tiny_step_size
241 : CASE (linesearch_method_gold)
242 44 : ALLOCATE (this%state_gold)
243 44 : this%state_gold%last_step_size = this%init_step_size
244 44 : this%state_gold%eps_step_size = this%eps_step_size
245 : CASE (linesearch_method_none)
246 : ! nothing todo
247 : CASE DEFAULT
248 526 : CPABORT("unknown method")
249 : END SELECT
250 :
251 526 : END SUBROUTINE linesearch_init_low
252 :
253 : ! **************************************************************************************************
254 : !> \brief Finzalize line search machinery
255 : !> \param this ...
256 : !> \author Ole Schuett
257 : ! **************************************************************************************************
258 526 : SUBROUTINE linesearch_finalize(this)
259 : TYPE(linesearch_type), INTENT(INOUT) :: this
260 :
261 928 : SELECT CASE (this%method)
262 : CASE (linesearch_method_adapt)
263 402 : DEALLOCATE (this%state_adapt)
264 : CASE (linesearch_method_2pnt)
265 16 : DEALLOCATE (this%state_2pnt)
266 : CASE (linesearch_method_3pnt)
267 60 : DEALLOCATE (this%state_3pnt)
268 : CASE (linesearch_method_gold)
269 44 : DEALLOCATE (this%state_gold)
270 : CASE (linesearch_method_none)
271 : ! nothing todo
272 : CASE DEFAULT
273 526 : CPABORT("unknown method")
274 : END SELECT
275 :
276 : !TODO: should finish printkey, but don't have the section here
277 526 : END SUBROUTINE linesearch_finalize
278 :
279 : ! **************************************************************************************************
280 : !> \brief Reset line search to initial state
281 : !> \param this ...
282 : !> \author Ole Schuett
283 : ! **************************************************************************************************
284 246 : SUBROUTINE linesearch_reset(this)
285 : TYPE(linesearch_type), INTENT(INOUT) :: this
286 :
287 246 : CALL linesearch_finalize(this)
288 246 : CALL linesearch_init_low(this)
289 246 : 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 9896 : 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 17322 : SELECT CASE (this%method)
306 : CASE (linesearch_method_adapt)
307 7426 : CALL linesearch_adapt(this%state_adapt, energy, step_size, is_done, this%iw, TRIM(this%label))
308 : CASE (linesearch_method_2pnt)
309 68 : CALL linesearch_2pnt(this%state_2pnt, energy, slope, step_size, is_done, this%iw, TRIM(this%label))
310 : CASE (linesearch_method_3pnt)
311 768 : CALL linesearch_3pnt(this%state_3pnt, energy, step_size, is_done, this%iw, TRIM(this%label))
312 : CASE (linesearch_method_gold)
313 1634 : CALL linesearch_gold(this%state_gold, energy, step_size, is_done, this%iw, TRIM(this%label))
314 : CASE (linesearch_method_none)
315 0 : step_size = this%init_step_size ! take steps of fixed length
316 0 : is_done = .TRUE.
317 : CASE DEFAULT
318 9896 : CPABORT("unknown method")
319 : END SELECT
320 :
321 9896 : this%step_size = step_size
322 9896 : this%starts = is_done
323 9896 : 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 68 : 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 68 : this%energies(this%count) = energy
347 68 : is_done = .FALSE.
348 :
349 34 : SELECT CASE (this%count)
350 : CASE (1)
351 34 : step_size = 2.0_dp*this%last_step_size
352 34 : this%scan_step = step_size
353 34 : this%count = 2
354 : CASE (2)
355 34 : c = this%energies(1)
356 34 : b = -slope
357 34 : x2 = this%scan_step
358 34 : a = (this%energies(2) - b*x2 - c)/(x2**2)
359 :
360 34 : IF (a < 0.0_dp) THEN
361 0 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| had to quench curvature"
362 : a = 1.0E-15_dp
363 : END IF
364 :
365 34 : step_size = -b/(2.0_dp*a)
366 34 : pred_energy = a*step_size**2 + b*step_size + c
367 :
368 34 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 2pnt suggested step_size: ", step_size
369 34 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 2pnt predicted energy", pred_energy
370 :
371 34 : IF (pred_energy > this%energies(1) .OR. pred_energy > this%energies(2)) THEN
372 0 : CPABORT(label//"LS| predicted energy not below test points")
373 : END IF
374 :
375 34 : IF (step_size > this%max_step_size) THEN
376 0 : step_size = this%max_step_size
377 0 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| limiting step_size to MAX_STEP_SIZE"
378 : END IF
379 :
380 34 : this%last_step_size = step_size
381 34 : this%count = 1
382 34 : is_done = .TRUE.
383 : CASE DEFAULT
384 68 : CPABORT("this should not happen")
385 : END SELECT
386 :
387 68 : 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 768 : 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 768 : this%energies(this%count) = energy
411 768 : is_done = .FALSE.
412 :
413 256 : SELECT CASE (this%count)
414 : CASE (1)
415 256 : step_size = (2.0_dp/3.0_dp)*this%last_step_size
416 256 : IF (step_size < this%tiny_step_size) THEN
417 0 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| initial step size too small, using TINY_STEP_SIZE"
418 0 : step_size = this%tiny_step_size
419 : END IF
420 256 : this%scan_steps(1) = 0.0_dp
421 256 : this%scan_steps(2) = step_size
422 256 : this%count = 2
423 : CASE (2)
424 256 : IF (this%energies(1) > this%energies(2)) THEN
425 216 : step_size = 2.0_dp*this%scan_steps(2)
426 : ELSE
427 40 : step_size = 0.5_dp*this%scan_steps(2)
428 : END IF
429 256 : this%scan_steps(3) = step_size
430 256 : this%count = 3
431 : CASE (3)
432 : ! fitting y = a*x^2 + b*x + c
433 256 : y1 = this%energies(1); y2 = this%energies(2); y3 = this%energies(3)
434 256 : x1 = this%scan_steps(1); x2 = this%scan_steps(2); x3 = this%scan_steps(3)
435 :
436 256 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt scan_steps: ", this%scan_steps
437 256 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt energies: ", this%energies
438 :
439 : ! Cramer's Rule
440 256 : denom = (x1 - x2)*(x1 - x3)*(x2 - x3)
441 256 : a = (x3*(y2 - y1) + x2*(y1 - y3) + x1*(y3 - y2))/denom
442 256 : b = (x3**2*(y1 - y2) + x2**2*(y3 - y1) + x1**2*(y2 - y3))/denom
443 256 : c = (x2*x3*(x2 - x3)*y1 + x3*x1*(x3 - x1)*y2 + x1*x2*(x1 - x2)*y3)/denom
444 :
445 256 : step_size = -b/(2.0_dp*a)
446 256 : pred_energy = a*step_size**2 + b*step_size + c
447 256 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt suggested step_size: ", step_size
448 256 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| 3pnt predicted energy", pred_energy
449 :
450 256 : IF (a < 0) THEN
451 0 : step_size = -2.0_dp*step_size
452 0 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| inverting step size"
453 : END IF
454 :
455 256 : IF (step_size < 0) THEN
456 0 : step_size = this%tiny_step_size
457 0 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| makeing a step of size TINY_STEP_SIZE"
458 : END IF
459 :
460 256 : IF (step_size > this%max_step_size) THEN
461 6 : step_size = this%max_step_size
462 6 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| limiting step_size to MAX_STEP_SIZE"
463 : END IF
464 :
465 256 : this%last_step_size = step_size
466 256 : this%count = 1
467 256 : is_done = .TRUE.
468 : CASE DEFAULT
469 768 : CPABORT("this should not happen")
470 : END SELECT
471 :
472 768 : 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 7426 : 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 7426 : is_done = .FALSE.
499 7426 : this%count = this%count + 1
500 :
501 7426 : IF (.NOT. this%have_left) THEN
502 1932 : this%left_x = 0.0_dp
503 1932 : this%left_e = energy
504 1932 : this%have_left = .TRUE.
505 1932 : step_size = this%last_step_size
506 :
507 5494 : ELSE IF (.NOT. (this%have_middle .OR. this%have_right)) THEN
508 1932 : IF (energy < this%left_e) THEN
509 1550 : this%middle_e = energy
510 1550 : this%middle_x = this%last_step_size
511 1550 : this%have_middle = .TRUE.
512 1550 : step_size = this%middle_x*grow_factor
513 : ELSE
514 382 : this%right_e = energy
515 382 : this%right_x = this%last_step_size
516 382 : this%have_right = .TRUE.
517 382 : step_size = this%right_x*shrink_factor
518 : END IF
519 :
520 3562 : ELSE IF (.NOT. this%have_right) THEN
521 2978 : IF (energy < this%middle_e) THEN
522 1430 : this%middle_e = energy
523 1430 : this%middle_x = this%last_step_size
524 1430 : step_size = this%middle_x*grow_factor
525 : ELSE
526 1548 : this%right_e = energy
527 1548 : this%right_x = this%last_step_size
528 1548 : this%have_right = .TRUE.
529 : END IF
530 :
531 584 : ELSE IF (.NOT. this%have_middle) THEN
532 584 : IF (energy > this%left_e) THEN
533 202 : this%right_e = energy
534 202 : this%right_x = this%last_step_size
535 202 : step_size = this%right_x*shrink_factor
536 : ELSE
537 382 : this%middle_e = energy
538 382 : this%middle_x = this%last_step_size
539 382 : this%have_middle = .TRUE.
540 : END IF
541 : END IF
542 :
543 7426 : IF (this%count > 3) THEN
544 1632 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| Need extra step"
545 : END IF
546 7426 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: ", this%have_left, this%have_middle, this%have_right
547 7426 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: scan_steps: ", this%left_x, this%middle_x, this%right_x
548 7426 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: energies: ", this%left_e, this%middle_e, this%right_e
549 :
550 7426 : IF (this%have_left .AND. this%have_middle .AND. this%have_right) THEN
551 : ! fitting y = a*x^2 + b*x + c
552 1930 : y1 = this%left_e; y2 = this%middle_e; y3 = this%right_e
553 1930 : x1 = this%left_x; x2 = this%middle_x; x3 = this%right_x
554 :
555 : ! Cramer's rule
556 1930 : denom = (x1 - x2)*(x1 - x3)*(x2 - x3)
557 1930 : a = (x3*(y2 - y1) + x2*(y1 - y3) + x1*(y3 - y2))/denom
558 1930 : b = (x3**2*(y1 - y2) + x2**2*(y3 - y1) + x1**2*(y2 - y3))/denom
559 1930 : c = (x2*x3*(x2 - x3)*y1 + x3*x1*(x3 - x1)*y2 + x1*x2*(x1 - x2)*y3)/denom
560 :
561 1930 : IF (ABS(a) /= 0.0_dp) THEN
562 1930 : step_size = -b/(2.0_dp*a)
563 : ELSE
564 0 : step_size = 0.0_dp
565 : END IF
566 :
567 1930 : CPASSERT(step_size >= 0.0_dp)
568 1930 : pred_energy = a*step_size**2 + b*step_size + c
569 1930 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: suggested step_size: ", step_size
570 1930 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| adapt: predicted energy", pred_energy
571 :
572 : ! reset
573 1930 : is_done = .TRUE.
574 1930 : this%count = 0
575 1930 : this%have_left = .FALSE.
576 1930 : this%have_middle = .FALSE.
577 1930 : this%have_right = .FALSE.
578 1930 : this%left_e = 0.0
579 1930 : this%middle_e = 0.0
580 1930 : this%right_e = 0.0
581 1930 : this%left_x = 0.0
582 1930 : this%middle_x = 0.0
583 1930 : this%right_x = 0.0
584 : END IF
585 :
586 7426 : this%last_step_size = step_size
587 7426 : 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 1634 : 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 1634 : is_done = .FALSE.
612 :
613 1634 : IF (this%gave_up) &
614 0 : CPABORT("had to give up, should not be called again")
615 :
616 1634 : IF (.NOT. this%have_left) THEN
617 154 : this%left_x = 0.0_dp
618 154 : this%left_e = energy
619 154 : this%have_left = .TRUE.
620 154 : step_size = this%last_step_size
621 :
622 1480 : ELSE IF (.NOT. (this%have_middle .OR. this%have_right)) THEN
623 154 : IF (energy < this%left_e) THEN
624 112 : this%middle_e = energy
625 112 : this%middle_x = this%scan_steps
626 112 : this%have_middle = .TRUE.
627 112 : step_size = this%middle_x*phi
628 : ELSE
629 42 : this%right_e = energy
630 42 : this%right_x = this%scan_steps
631 42 : this%have_right = .TRUE.
632 42 : step_size = this%right_x/phi
633 : END IF
634 :
635 1326 : ELSE IF (.NOT. this%have_right) THEN
636 418 : IF (energy < this%middle_e) THEN
637 306 : this%middle_e = energy
638 306 : this%middle_x = this%scan_steps
639 306 : step_size = this%middle_x*phi
640 : ELSE
641 112 : this%right_e = energy
642 112 : this%right_x = this%scan_steps
643 112 : this%have_right = .TRUE.
644 : END IF
645 :
646 908 : ELSE IF (.NOT. this%have_middle) THEN
647 112 : IF (energy > this%left_e) THEN
648 70 : this%right_e = energy
649 70 : this%right_x = this%scan_steps
650 70 : step_size = this%right_x/phi
651 : ELSE
652 42 : this%middle_e = energy
653 42 : this%middle_x = this%scan_steps
654 42 : this%have_middle = .TRUE.
655 : END IF
656 :
657 : ELSE !up and running
658 796 : a = this%middle_x - this%left_x
659 796 : b = this%right_x - this%middle_x
660 796 : IF (energy < this%middle_e) THEN
661 234 : IF (a < b) THEN
662 90 : this%left_e = this%middle_e
663 90 : this%left_x = this%middle_x
664 : ELSE
665 144 : this%right_e = this%middle_e
666 144 : this%right_x = this%middle_x
667 : END IF
668 234 : this%middle_e = energy
669 234 : this%middle_x = this%scan_steps
670 : ELSE
671 562 : IF (a < b) THEN
672 234 : this%right_e = energy
673 234 : this%right_x = this%scan_steps
674 : ELSE
675 328 : this%left_e = energy
676 328 : this%left_x = this%scan_steps
677 : END IF
678 : END IF
679 : END IF
680 :
681 1634 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%have_left, this%have_middle, this%have_right
682 1634 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%left_x, this%middle_x, this%right_x
683 1634 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold: ", this%left_e, this%middle_e, this%right_e
684 :
685 1634 : IF (this%have_left .AND. this%have_middle .AND. this%have_right) THEN
686 950 : a = this%middle_x - this%left_x
687 950 : b = this%right_x - this%middle_x
688 950 : IF (ABS(MIN(a, b)*phi - MAX(a, b)) > 1.0E-10) &
689 0 : CPABORT("golden-ratio gone")
690 :
691 950 : IF (a < b) THEN
692 418 : step_size = this%middle_x + a/phi
693 : ELSE
694 532 : step_size = this%middle_x - b/phi
695 : END IF
696 :
697 950 : d = ABS(this%right_x - this%left_x)/(ABS(this%middle_x) + ABS(step_size))
698 950 : IF (d < this%eps_step_size) THEN
699 154 : step_size = this%middle_x
700 154 : this%last_step_size = step_size
701 154 : is_done = .TRUE.
702 :
703 154 : IF (unit_nr > 0) WRITE (unit_nr, *) label, "LS| gold done, step-size: ", step_size
704 :
705 154 : this%have_left = .FALSE.
706 154 : this%have_middle = .FALSE.
707 154 : this%have_right = .FALSE.
708 154 : this%left_e = 0.0
709 154 : this%middle_e = 0.0
710 154 : this%right_e = 0.0
711 154 : this%left_x = 0.0
712 154 : this%middle_x = 0.0
713 154 : this%right_x = 0.0
714 : END IF
715 :
716 : END IF
717 :
718 1634 : IF (step_size < 1E-10) CPABORT("linesearch failed / done")
719 :
720 1634 : this%scan_steps = step_size
721 1634 : END SUBROUTINE linesearch_gold
722 :
723 0 : END MODULE linesearch
|
