Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines that optimize a functional using the limited memory bfgs
10 : !> quasi-newton method.
11 : !> The process set up so that a master runs the real optimizer and the
12 : !> others help then to calculate the objective function.
13 : !> The arguments for the objective function are physically present in
14 : !> every processor (nedeed in the actual implementation of pao).
15 : !> In the future tha arguments themselves could be distributed.
16 : !> \par History
17 : !> 09.2003 globenv->para_env, retain/release, better parallel behaviour
18 : !> 01.2020 Space Group Symmetry introduced by Pierre-André Cazade [pcazade]
19 : !> \author Fawzi Mohamed
20 : !> @version 2.2002
21 : ! **************************************************************************************************
22 : MODULE cp_lbfgs_optimizer_gopt
23 : USE cp_lbfgs, ONLY: setulb
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_type,&
26 : cp_to_string
27 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
28 : cp_print_key_unit_nr
29 : USE cp_subsys_types, ONLY: cp_subsys_type
30 : USE force_env_types, ONLY: force_env_get,&
31 : force_env_type
32 : USE gopt_f_methods, ONLY: cp_eval_at,&
33 : gopt_f_io
34 : USE gopt_f_types, ONLY: gopt_f_release,&
35 : gopt_f_retain,&
36 : gopt_f_type
37 : USE gopt_param_types, ONLY: gopt_param_type
38 : USE input_section_types, ONLY: section_vals_type
39 : USE kinds, ONLY: dp
40 : USE machine, ONLY: m_walltime
41 : USE message_passing, ONLY: mp_para_env_release,&
42 : mp_para_env_type
43 : USE space_groups, ONLY: spgr_apply_rotations_coord,&
44 : spgr_apply_rotations_force
45 : USE space_groups_types, ONLY: spgr_type
46 : #include "../base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 : PRIVATE
50 :
51 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_lbfgs_optimizer_gopt'
53 :
54 : ! types
55 : PUBLIC :: cp_lbfgs_opt_gopt_type
56 :
57 : ! core methods
58 :
59 : ! special methos
60 :
61 : ! underlying functions
62 : PUBLIC :: cp_opt_gopt_create, cp_opt_gopt_release, &
63 : cp_opt_gopt_next, &
64 : cp_opt_gopt_stop
65 :
66 : ! **************************************************************************************************
67 : !> \brief info for the optimizer (see the description of this module)
68 : !> \param task the actual task of the optimizer (in the master it is up to
69 : !> date, in case of error also the minions one get updated.
70 : !> \param csave internal character string used by the lbfgs optimizer,
71 : !> meaningful only in the master
72 : !> \param lsave logical array used by the lbfgs optimizer, updated only
73 : !> in the master
74 : !> On exit with task = 'NEW_X', the following information is
75 : !> available:
76 : !> lsave(1) = .true. the initial x did not satisfy the bounds;
77 : !> lsave(2) = .true. the problem contains bounds;
78 : !> lsave(3) = .true. each variable has upper and lower bounds.
79 : !> \param ref_count reference count (see doc/ReferenceCounting.html)
80 : !> \param m the dimension of the subspace used to approximate the second
81 : !> derivative
82 : !> \param print_every every how many iterations output should be written.
83 : !> if 0 only at end, if print_every<0 never
84 : !> \param master the pid of the master processor
85 : !> \param max_f_per_iter the maximum number of function evaluations per
86 : !> iteration
87 : !> \param status 0: just initialized, 1: f g calculation,
88 : !> 2: begin new iteration, 3: ended iteration,
89 : !> 4: normal (converged) exit, 5: abnormal (error) exit,
90 : !> 6: daellocated
91 : !> \param n_iter the actual iteration number
92 : !> \param kind_of_bound an array with 0 (no bound), 1 (lower bound),
93 : !> 2 (both bounds), 3 (upper bound), to describe the bounds
94 : !> of every variable
95 : !> \param i_work_array an integer workarray of dimension 3*n, present only
96 : !> in the master
97 : !> \param isave is an INTEGER working array of dimension 44.
98 : !> On exit with task = 'NEW_X', it contains information that
99 : !> the user may want to access:
100 : !> \param isave (30) = the current iteration number;
101 : !> \param isave (34) = the total number of function and gradient
102 : !> evaluations;
103 : !> \param isave (36) = the number of function value or gradient
104 : !> evaluations in the current iteration;
105 : !> \param isave (38) = the number of free variables in the current
106 : !> iteration;
107 : !> \param isave (39) = the number of active constraints at the current
108 : !> iteration;
109 : !> \param f the actual best value of the object function
110 : !> \param wanted_relative_f_delta the wanted relative error on f
111 : !> (to be multiplied by epsilon), 0.0 -> no check
112 : !> \param wanted_projected_gradient the wanted error on the projected
113 : !> gradient (hessian times the gradient), 0.0 -> no check
114 : !> \param last_f the value of f in the last iteration
115 : !> \param projected_gradient the value of the sup norm of the projected
116 : !> gradient
117 : !> \param x the actual evaluation point (best one if converged or stopped)
118 : !> \param lower_bound the lower bounds
119 : !> \param upper_bound the upper bounds
120 : !> \param gradient the actual gradient
121 : !> \param dsave info date for lbfgs (master only)
122 : !> \param work_array a work array for lbfgs (master only)
123 : !> \param para_env the parallel environment for this optimizer
124 : !> \param obj_funct the objective function to be optimized
125 : !> \par History
126 : !> none
127 : !> \author Fawzi Mohamed
128 : !> @version 2.2002
129 : ! **************************************************************************************************
130 : TYPE cp_lbfgs_opt_gopt_type
131 : CHARACTER(len=60) :: task = ""
132 : CHARACTER(len=60) :: csave = ""
133 : LOGICAL :: lsave(4) = .FALSE.
134 : INTEGER :: m = 0, print_every = 0, master = 0, max_f_per_iter = 0, status = 0, n_iter = 0
135 : INTEGER, DIMENSION(:), POINTER :: kind_of_bound => NULL(), i_work_array => NULL(), isave => NULL()
136 : REAL(kind=dp) :: f = 0.0_dp, wanted_relative_f_delta = 0.0_dp, wanted_projected_gradient = 0.0_dp, &
137 : last_f = 0.0_dp, projected_gradient = 0.0_dp, eold = 0.0_dp, emin = 0.0_dp, trust_radius = 0.0_dp
138 : REAL(kind=dp), DIMENSION(:), POINTER :: x => NULL(), lower_bound => NULL(), upper_bound => NULL(), &
139 : gradient => NULL(), dsave => NULL(), work_array => NULL()
140 : TYPE(mp_para_env_type), POINTER :: para_env => NULL()
141 : TYPE(gopt_f_type), POINTER :: obj_funct => NULL()
142 : END TYPE cp_lbfgs_opt_gopt_type
143 :
144 : CONTAINS
145 :
146 : ! **************************************************************************************************
147 : !> \brief calls the L-BFGS optimizer
148 : !> \param optimizer the optimizer state
149 : !> \param spgr optional space group information
150 : !> \param iwunit optional output unit
151 : ! **************************************************************************************************
152 3192 : SUBROUTINE cp_opt_gopt_setulb(optimizer, spgr, iwunit)
153 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
154 : TYPE(spgr_type), OPTIONAL, POINTER :: spgr
155 : INTEGER, OPTIONAL :: iwunit
156 :
157 : CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
158 : optimizer%lower_bound, optimizer%upper_bound, &
159 : optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
160 : optimizer%wanted_relative_f_delta, &
161 : optimizer%wanted_projected_gradient, optimizer%work_array, &
162 : optimizer%i_work_array, optimizer%task, optimizer%print_every, &
163 : optimizer%csave, optimizer%lsave, optimizer%isave, &
164 3192 : optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=iwunit)
165 :
166 3192 : END SUBROUTINE cp_opt_gopt_setulb
167 :
168 : ! **************************************************************************************************
169 : !> \brief initializes the optimizer
170 : !> \param optimizer ...
171 : !> \param para_env ...
172 : !> \param obj_funct ...
173 : !> \param x0 ...
174 : !> \param m ...
175 : !> \param print_every ...
176 : !> \param wanted_relative_f_delta ...
177 : !> \param wanted_projected_gradient ...
178 : !> \param lower_bound ...
179 : !> \param upper_bound ...
180 : !> \param kind_of_bound ...
181 : !> \param master ...
182 : !> \param max_f_per_iter ...
183 : !> \param trust_radius ...
184 : !> \par History
185 : !> 02.2002 created [fawzi]
186 : !> 09.2003 refactored (retain/release,para_env) [fawzi]
187 : !> \author Fawzi Mohamed
188 : !> \note
189 : !> redirects the lbfgs output the the default unit
190 : ! **************************************************************************************************
191 540 : SUBROUTINE cp_opt_gopt_create(optimizer, para_env, obj_funct, x0, m, print_every, &
192 0 : wanted_relative_f_delta, wanted_projected_gradient, lower_bound, upper_bound, &
193 0 : kind_of_bound, master, max_f_per_iter, trust_radius)
194 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(OUT) :: optimizer
195 : TYPE(mp_para_env_type), POINTER :: para_env
196 : TYPE(gopt_f_type), POINTER :: obj_funct
197 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: x0
198 : INTEGER, INTENT(in), OPTIONAL :: m, print_every
199 : REAL(kind=dp), INTENT(in), OPTIONAL :: wanted_relative_f_delta, &
200 : wanted_projected_gradient
201 : REAL(kind=dp), DIMENSION(SIZE(x0)), INTENT(in), &
202 : OPTIONAL :: lower_bound, upper_bound
203 : INTEGER, DIMENSION(SIZE(x0)), INTENT(in), OPTIONAL :: kind_of_bound
204 : INTEGER, INTENT(in), OPTIONAL :: master, max_f_per_iter
205 : REAL(kind=dp), INTENT(in), OPTIONAL :: trust_radius
206 :
207 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_create'
208 :
209 : INTEGER :: handle, lenwa, n
210 :
211 90 : CALL timeset(routineN, handle)
212 :
213 : NULLIFY (optimizer%kind_of_bound, &
214 90 : optimizer%i_work_array, &
215 90 : optimizer%isave, &
216 90 : optimizer%x, &
217 90 : optimizer%lower_bound, &
218 90 : optimizer%upper_bound, &
219 90 : optimizer%gradient, &
220 90 : optimizer%dsave, &
221 90 : optimizer%work_array, &
222 : optimizer%para_env, &
223 90 : optimizer%obj_funct)
224 90 : n = SIZE(x0)
225 90 : optimizer%m = 4
226 90 : IF (PRESENT(m)) optimizer%m = m
227 90 : optimizer%master = para_env%source
228 90 : optimizer%para_env => para_env
229 90 : CALL para_env%retain()
230 90 : optimizer%obj_funct => obj_funct
231 90 : CALL gopt_f_retain(obj_funct)
232 90 : optimizer%max_f_per_iter = 20
233 90 : IF (PRESENT(max_f_per_iter)) optimizer%max_f_per_iter = max_f_per_iter
234 90 : optimizer%print_every = -1
235 90 : optimizer%n_iter = 0
236 90 : optimizer%f = -1.0_dp
237 90 : optimizer%last_f = -1.0_dp
238 90 : optimizer%projected_gradient = -1.0_dp
239 90 : IF (PRESENT(print_every)) optimizer%print_every = print_every
240 90 : IF (PRESENT(master)) optimizer%master = master
241 90 : IF (optimizer%master == optimizer%para_env%mepos) THEN
242 : !MK This has to be adapted for a new L-BFGS version possibly
243 45 : lenwa = 2*optimizer%m*n + 5*n + 11*optimizer%m*optimizer%m + 8*optimizer%m
244 : ALLOCATE (optimizer%kind_of_bound(n), optimizer%i_work_array(3*n), &
245 225 : optimizer%isave(44))
246 : ALLOCATE (optimizer%x(n), optimizer%lower_bound(n), &
247 : optimizer%upper_bound(n), optimizer%gradient(n), &
248 360 : optimizer%dsave(29), optimizer%work_array(lenwa))
249 28632 : optimizer%x = x0
250 45 : optimizer%task = 'START'
251 85806 : optimizer%i_work_array = 0
252 2025 : optimizer%isave = 0
253 28632 : optimizer%lower_bound = 0.0_dp
254 28632 : optimizer%upper_bound = 0.0_dp
255 28632 : optimizer%gradient = 0.0_dp
256 1350 : optimizer%dsave = 0.0_dp
257 544950 : optimizer%work_array = 0.0_dp
258 45 : IF (PRESENT(wanted_relative_f_delta)) THEN
259 45 : optimizer%wanted_relative_f_delta = wanted_relative_f_delta
260 : END IF
261 45 : IF (PRESENT(wanted_projected_gradient)) THEN
262 45 : optimizer%wanted_projected_gradient = wanted_projected_gradient
263 : END IF
264 28632 : optimizer%kind_of_bound = 0
265 45 : IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound
266 45 : IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound
267 45 : IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound
268 45 : IF (PRESENT(trust_radius)) optimizer%trust_radius = trust_radius
269 :
270 45 : CALL cp_opt_gopt_setulb(optimizer)
271 : ELSE
272 : NULLIFY ( &
273 45 : optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, &
274 45 : optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, &
275 45 : optimizer%dsave, optimizer%work_array)
276 135 : ALLOCATE (optimizer%x(n))
277 28632 : optimizer%x(:) = 0.0_dp
278 90 : ALLOCATE (optimizer%gradient(n))
279 28632 : optimizer%gradient(:) = 0.0_dp
280 : END IF
281 114438 : CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
282 90 : optimizer%status = 0
283 :
284 90 : CALL timestop(handle)
285 :
286 90 : END SUBROUTINE cp_opt_gopt_create
287 :
288 : ! **************************************************************************************************
289 : !> \brief releases the optimizer (see doc/ReferenceCounting.html)
290 : !> \param optimizer the object that should be released
291 : !> \par History
292 : !> 02.2002 created [fawzi]
293 : !> 09.2003 dealloc_ref->release [fawzi]
294 : !> \author Fawzi Mohamed
295 : ! **************************************************************************************************
296 90 : SUBROUTINE cp_opt_gopt_release(optimizer)
297 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
298 :
299 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_release'
300 :
301 : INTEGER :: handle
302 :
303 90 : CALL timeset(routineN, handle)
304 :
305 90 : IF (ASSOCIATED(optimizer%kind_of_bound)) THEN
306 45 : DEALLOCATE (optimizer%kind_of_bound)
307 : END IF
308 90 : IF (ASSOCIATED(optimizer%i_work_array)) THEN
309 45 : DEALLOCATE (optimizer%i_work_array)
310 : END IF
311 90 : IF (ASSOCIATED(optimizer%isave)) THEN
312 45 : DEALLOCATE (optimizer%isave)
313 : END IF
314 90 : IF (ASSOCIATED(optimizer%x)) THEN
315 90 : DEALLOCATE (optimizer%x)
316 : END IF
317 90 : IF (ASSOCIATED(optimizer%lower_bound)) THEN
318 45 : DEALLOCATE (optimizer%lower_bound)
319 : END IF
320 90 : IF (ASSOCIATED(optimizer%upper_bound)) THEN
321 45 : DEALLOCATE (optimizer%upper_bound)
322 : END IF
323 90 : IF (ASSOCIATED(optimizer%gradient)) THEN
324 90 : DEALLOCATE (optimizer%gradient)
325 : END IF
326 90 : IF (ASSOCIATED(optimizer%dsave)) THEN
327 45 : DEALLOCATE (optimizer%dsave)
328 : END IF
329 90 : IF (ASSOCIATED(optimizer%work_array)) THEN
330 45 : DEALLOCATE (optimizer%work_array)
331 : END IF
332 90 : CALL mp_para_env_release(optimizer%para_env)
333 90 : CALL gopt_f_release(optimizer%obj_funct)
334 :
335 90 : CALL timestop(handle)
336 90 : END SUBROUTINE cp_opt_gopt_release
337 :
338 : ! **************************************************************************************************
339 : !> \brief takes different valuse from the optimizer
340 : !> \param optimizer ...
341 : !> \param para_env ...
342 : !> \param obj_funct ...
343 : !> \param m ...
344 : !> \param print_every ...
345 : !> \param wanted_relative_f_delta ...
346 : !> \param wanted_projected_gradient ...
347 : !> \param x ...
348 : !> \param lower_bound ...
349 : !> \param upper_bound ...
350 : !> \param kind_of_bound ...
351 : !> \param master ...
352 : !> \param actual_projected_gradient ...
353 : !> \param n_var ...
354 : !> \param n_iter ...
355 : !> \param status ...
356 : !> \param max_f_per_iter ...
357 : !> \param at_end ...
358 : !> \param is_master ...
359 : !> \param last_f ...
360 : !> \param f ...
361 : !> \par History
362 : !> none
363 : !> \author Fawzi Mohamed
364 : !> @version 2.2002
365 : ! **************************************************************************************************
366 0 : SUBROUTINE cp_opt_gopt_get(optimizer, para_env, &
367 : obj_funct, m, print_every, &
368 : wanted_relative_f_delta, wanted_projected_gradient, &
369 : x, lower_bound, upper_bound, kind_of_bound, master, &
370 : actual_projected_gradient, &
371 : n_var, n_iter, status, max_f_per_iter, at_end, &
372 : is_master, last_f, f)
373 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
374 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
375 : TYPE(gopt_f_type), OPTIONAL, POINTER :: obj_funct
376 : INTEGER, INTENT(out), OPTIONAL :: m, print_every
377 : REAL(kind=dp), INTENT(out), OPTIONAL :: wanted_relative_f_delta, &
378 : wanted_projected_gradient
379 : REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: x, lower_bound, upper_bound
380 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: kind_of_bound
381 : INTEGER, INTENT(out), OPTIONAL :: master
382 : REAL(kind=dp), INTENT(out), OPTIONAL :: actual_projected_gradient
383 : INTEGER, INTENT(out), OPTIONAL :: n_var, n_iter, status, max_f_per_iter
384 : LOGICAL, INTENT(out), OPTIONAL :: at_end, is_master
385 : REAL(kind=dp), INTENT(out), OPTIONAL :: last_f, f
386 :
387 0 : IF (PRESENT(is_master)) is_master = optimizer%master == optimizer%para_env%mepos
388 0 : IF (PRESENT(master)) master = optimizer%master
389 0 : IF (PRESENT(status)) status = optimizer%status
390 0 : IF (PRESENT(para_env)) para_env => optimizer%para_env
391 0 : IF (PRESENT(obj_funct)) obj_funct = optimizer%obj_funct
392 0 : IF (PRESENT(m)) m = optimizer%m
393 0 : IF (PRESENT(max_f_per_iter)) max_f_per_iter = optimizer%max_f_per_iter
394 0 : IF (PRESENT(wanted_projected_gradient)) THEN
395 0 : wanted_projected_gradient = optimizer%wanted_projected_gradient
396 : END IF
397 0 : IF (PRESENT(wanted_relative_f_delta)) THEN
398 0 : wanted_relative_f_delta = optimizer%wanted_relative_f_delta
399 : END IF
400 0 : IF (PRESENT(print_every)) print_every = optimizer%print_every
401 0 : IF (PRESENT(x)) x => optimizer%x
402 0 : IF (PRESENT(n_var)) n_var = SIZE(x)
403 0 : IF (PRESENT(lower_bound)) lower_bound => optimizer%lower_bound
404 0 : IF (PRESENT(upper_bound)) upper_bound => optimizer%upper_bound
405 0 : IF (PRESENT(kind_of_bound)) kind_of_bound => optimizer%kind_of_bound
406 0 : IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
407 0 : IF (PRESENT(last_f)) last_f = optimizer%last_f
408 0 : IF (PRESENT(f)) f = optimizer%f
409 0 : IF (PRESENT(at_end)) at_end = optimizer%status > 3
410 0 : IF (PRESENT(actual_projected_gradient)) THEN
411 0 : actual_projected_gradient = optimizer%projected_gradient
412 : END IF
413 0 : IF (optimizer%master == optimizer%para_env%mepos) THEN
414 0 : IF (optimizer%isave(30) > 1 .AND. (optimizer%task(1:5) == "NEW_X" .OR. &
415 : optimizer%task(1:4) == "STOP" .AND. optimizer%task(7:9) == "CPU")) THEN
416 : ! nr iterations >1 .and. dsave contains the wanted data
417 0 : IF (PRESENT(last_f)) last_f = optimizer%dsave(2)
418 0 : IF (PRESENT(actual_projected_gradient)) THEN
419 0 : actual_projected_gradient = optimizer%dsave(13)
420 : END IF
421 : ELSE
422 0 : CPASSERT(.NOT. PRESENT(last_f))
423 0 : CPASSERT(.NOT. PRESENT(actual_projected_gradient))
424 : END IF
425 0 : ELSE IF (PRESENT(lower_bound) .OR. PRESENT(upper_bound) .OR. PRESENT(kind_of_bound)) THEN
426 0 : CPWARN("asked undefined types")
427 : END IF
428 :
429 0 : END SUBROUTINE cp_opt_gopt_get
430 :
431 : ! **************************************************************************************************
432 : !> \brief does one optimization step
433 : !> \param optimizer ...
434 : !> \param n_iter ...
435 : !> \param f ...
436 : !> \param last_f ...
437 : !> \param projected_gradient ...
438 : !> \param converged ...
439 : !> \param geo_section ...
440 : !> \param force_env ...
441 : !> \param gopt_param ...
442 : !> \param spgr ...
443 : !> \par History
444 : !> 01.2020 modified [pcazade]
445 : !> \author Fawzi Mohamed
446 : !> @version 2.2002
447 : !> \note
448 : !> use directly mainlb in place of setulb ??
449 : ! **************************************************************************************************
450 2982 : SUBROUTINE cp_opt_gopt_step(optimizer, n_iter, f, last_f, &
451 : projected_gradient, converged, geo_section, force_env, &
452 : gopt_param, spgr)
453 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
454 : INTEGER, INTENT(out), OPTIONAL :: n_iter
455 : REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
456 : LOGICAL, INTENT(out), OPTIONAL :: converged
457 : TYPE(section_vals_type), POINTER :: geo_section
458 : TYPE(force_env_type), POINTER :: force_env
459 : TYPE(gopt_param_type), POINTER :: gopt_param
460 : TYPE(spgr_type), OPTIONAL, POINTER :: spgr
461 :
462 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_step'
463 :
464 : CHARACTER(LEN=5) :: wildcard
465 : INTEGER :: dataunit, handle, its
466 : LOGICAL :: conv, is_master, justEntred, &
467 : keep_space_group
468 : REAL(KIND=dp) :: t_diff, t_now, t_old
469 2982 : REAL(KIND=dp), DIMENSION(:), POINTER :: xold
470 : TYPE(cp_logger_type), POINTER :: logger
471 : TYPE(cp_subsys_type), POINTER :: subsys
472 :
473 2982 : NULLIFY (logger, xold)
474 5964 : logger => cp_get_default_logger()
475 2982 : CALL timeset(routineN, handle)
476 2982 : justEntred = .TRUE.
477 2982 : is_master = optimizer%master == optimizer%para_env%mepos
478 2982 : IF (PRESENT(converged)) converged = optimizer%status == 4
479 8946 : ALLOCATE (xold(SIZE(optimizer%x)))
480 :
481 : ! collecting subsys
482 2982 : CALL force_env_get(force_env, subsys=subsys)
483 :
484 2982 : keep_space_group = .FALSE.
485 2982 : IF (PRESENT(spgr)) THEN
486 2982 : IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
487 : END IF
488 :
489 : ! applies rotation matrices to coordinates
490 2982 : IF (keep_space_group) THEN
491 2 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
492 : END IF
493 :
494 1739322 : xold = optimizer%x
495 2982 : t_old = m_walltime()
496 :
497 2982 : IF (optimizer%status >= 4) THEN
498 0 : CPWARN("status>=4, trying to restart")
499 0 : optimizer%status = 0
500 : dataunit = cp_print_key_unit_nr(logger, geo_section, &
501 0 : "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
502 0 : IF (is_master) THEN
503 0 : optimizer%task = 'START'
504 0 : CALL cp_opt_gopt_setulb(optimizer, spgr=spgr, iwunit=dataunit)
505 : END IF
506 : CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
507 0 : "PRINT%PROGRAM_RUN_INFO")
508 : END IF
509 :
510 : DO
511 : dataunit = cp_print_key_unit_nr(logger, geo_section, &
512 9276 : "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
513 9276 : ifMaster: IF (is_master) THEN
514 4638 : IF (optimizer%task(1:7) == 'RESTART') THEN
515 : ! restart the optimizer
516 0 : optimizer%status = 0
517 0 : optimizer%task = 'START'
518 : ! applies rotation matrices to coordinates and forces
519 0 : IF (keep_space_group) THEN
520 0 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
521 0 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
522 : END IF
523 0 : CALL cp_opt_gopt_setulb(optimizer, spgr=spgr, iwunit=dataunit)
524 0 : IF (keep_space_group) THEN
525 0 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
526 0 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
527 : END IF
528 : END IF
529 4638 : IF (optimizer%task(1:2) == 'FG') THEN
530 1700 : IF (optimizer%isave(36) > optimizer%max_f_per_iter) THEN
531 0 : optimizer%task = 'STOP: CPU, hit max f eval in iter'
532 0 : optimizer%status = 5 ! anormal exit
533 0 : CALL cp_opt_gopt_setulb(optimizer, spgr=spgr, iwunit=dataunit)
534 : ELSE
535 1700 : optimizer%status = 1
536 : END IF
537 2938 : ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN
538 2937 : IF (justEntred) THEN
539 1447 : optimizer%status = 2
540 : ! applies rotation matrices to coordinates and forces
541 1447 : IF (keep_space_group) THEN
542 0 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
543 0 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
544 : END IF
545 1447 : CALL cp_opt_gopt_setulb(optimizer, spgr=spgr, iwunit=dataunit)
546 1447 : IF (keep_space_group) THEN
547 0 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
548 0 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
549 : END IF
550 : ELSE
551 : ! applies rotation matrices to coordinates and forces
552 1490 : IF (keep_space_group) THEN
553 1 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
554 1 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
555 : END IF
556 1490 : optimizer%status = 3
557 : END IF
558 1 : ELSE IF (optimizer%task(1:4) == 'CONV') THEN
559 1 : optimizer%status = 4
560 0 : ELSE IF (optimizer%task(1:4) == 'STOP') THEN
561 0 : optimizer%status = 5
562 0 : CPWARN("task became stop in an unknown way")
563 0 : ELSE IF (optimizer%task(1:5) == 'ERROR') THEN
564 0 : optimizer%status = 5
565 : ELSE
566 0 : CPWARN("unknown task '"//optimizer%task//"'")
567 : END IF
568 : END IF ifMaster
569 : CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
570 9276 : "PRINT%PROGRAM_RUN_INFO")
571 9276 : CALL optimizer%para_env%bcast(optimizer%status, optimizer%master)
572 : ! Dump info
573 9276 : IF (optimizer%status == 3) THEN
574 2980 : its = 0
575 2980 : IF (is_master) THEN
576 : ! Iteration level is taken into account in the optimizer external loop
577 1490 : its = optimizer%isave(30)
578 : END IF
579 : END IF
580 : !
581 3400 : SELECT CASE (optimizer%status)
582 : CASE (1)
583 : !op=1 evaluate f and g
584 : CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, &
585 : f=optimizer%f, &
586 : gradient=optimizer%gradient, &
587 3400 : master=optimizer%master, para_env=optimizer%para_env)
588 : ! do not use keywords?
589 : dataunit = cp_print_key_unit_nr(logger, geo_section, &
590 3400 : "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
591 3400 : IF (is_master) THEN
592 : ! applies rotation matrices to coordinates and forces
593 1700 : IF (keep_space_group) THEN
594 2 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
595 2 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
596 : END IF
597 1700 : CALL cp_opt_gopt_setulb(optimizer, spgr=spgr, iwunit=dataunit)
598 1700 : IF (keep_space_group) THEN
599 2 : CALL spgr_apply_rotations_coord(spgr, optimizer%x)
600 2 : CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
601 : END IF
602 : END IF
603 : CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
604 3400 : "PRINT%PROGRAM_RUN_INFO")
605 4131232 : CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
606 : CASE (2)
607 : !op=2 begin new iter
608 3361274 : CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
609 2894 : t_old = m_walltime()
610 : CASE (3)
611 : !op=3 ended iter
612 2980 : wildcard = "LBFGS"
613 : dataunit = cp_print_key_unit_nr(logger, geo_section, &
614 2980 : "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
615 2980 : IF (is_master) its = optimizer%isave(30)
616 2980 : CALL optimizer%para_env%bcast(its, optimizer%master)
617 :
618 : ! Some IO and Convergence check
619 2980 : t_now = m_walltime()
620 2980 : t_diff = t_now - t_old
621 2980 : t_old = t_now
622 : CALL gopt_f_io(optimizer%obj_funct, force_env, force_env%root_section, &
623 : its, optimizer%f, dataunit, optimizer%eold, optimizer%emin, wildcard, gopt_param, &
624 1739284 : SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff)
625 2980 : CALL optimizer%para_env%bcast(conv, optimizer%master)
626 : CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
627 2980 : "PRINT%PROGRAM_RUN_INFO")
628 2980 : optimizer%eold = optimizer%f
629 2980 : optimizer%emin = MIN(optimizer%emin, optimizer%eold)
630 1739284 : xold = optimizer%x
631 2980 : IF (PRESENT(converged)) converged = conv
632 2 : EXIT
633 : CASE (4)
634 : !op=4 (convergence - normal exit)
635 : ! Specific L-BFGS convergence criteria.. overrides the convergence criteria on
636 : ! stepsize and gradients
637 : dataunit = cp_print_key_unit_nr(logger, geo_section, &
638 2 : "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
639 2 : IF (dataunit > 0) THEN
640 1 : WRITE (dataunit, '(T2,A)') ""
641 1 : WRITE (dataunit, '(T2,A)') "************************************************"
642 1 : WRITE (dataunit, '(T2,A)') "* Specific L-BFGS convergence criteria *"
643 1 : WRITE (dataunit, '(T2,A)') "* WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR *"
644 1 : WRITE (dataunit, '(T2,A)') "* satisfied .... run CONVERGED! *"
645 1 : WRITE (dataunit, '(T2,A)') "* * * * *"
646 1 : WRITE (dataunit, '(T2,A)') "* General convergence criteria on stepsize and *"
647 1 : WRITE (dataunit, '(T2,A)') "* gradients may or may not have been satisfied *"
648 1 : WRITE (dataunit, '(T2,A)') "* yet; if unsatisfactory, try tightening the *"
649 1 : WRITE (dataunit, '(T2,A)') "* L-BFGS convergence criteria and restart run. *"
650 1 : WRITE (dataunit, '(T2,A)') "************************************************"
651 1 : WRITE (dataunit, '(T2,A)') ""
652 : END IF
653 : CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
654 2 : "PRINT%PROGRAM_RUN_INFO")
655 2 : IF (PRESENT(converged)) converged = .TRUE.
656 0 : EXIT
657 : CASE (5)
658 : ! Restore the last accepted point; a failed FG request may leave x at a trial point.
659 0 : CALL optimizer%para_env%bcast(optimizer%task, optimizer%master)
660 0 : optimizer%x = xold
661 : CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, &
662 : f=optimizer%f, gradient=optimizer%gradient, &
663 0 : master=optimizer%master, para_env=optimizer%para_env)
664 0 : IF (PRESENT(converged)) converged = .FALSE.
665 0 : EXIT
666 : CASE (6)
667 : ! deallocated
668 0 : CPABORT("step on a deallocated opt structure ")
669 : CASE default
670 : CALL cp_abort(__LOCATION__, &
671 0 : "unknown status "//cp_to_string(optimizer%status))
672 0 : optimizer%status = 5
673 9276 : EXIT
674 : END SELECT
675 6294 : IF (optimizer%status == 1 .AND. justEntred) THEN
676 88 : optimizer%eold = optimizer%f
677 88 : optimizer%emin = optimizer%eold
678 : END IF
679 : justEntred = .FALSE.
680 : END DO
681 :
682 3475662 : CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
683 : CALL cp_opt_gopt_bcast_res(optimizer, &
684 : n_iter=optimizer%n_iter, &
685 : f=optimizer%f, last_f=optimizer%last_f, &
686 2982 : projected_gradient=optimizer%projected_gradient)
687 :
688 2982 : DEALLOCATE (xold)
689 2982 : IF (PRESENT(f)) f = optimizer%f
690 2982 : IF (PRESENT(last_f)) last_f = optimizer%last_f
691 2982 : IF (PRESENT(projected_gradient)) projected_gradient = optimizer%projected_gradient
692 2982 : IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
693 2982 : CALL timestop(handle)
694 :
695 2982 : END SUBROUTINE cp_opt_gopt_step
696 :
697 : ! **************************************************************************************************
698 : !> \brief returns the results (and broadcasts them)
699 : !> \param optimizer the optimizer object the info is taken from
700 : !> \param n_iter the number of iterations
701 : !> \param f the actual value of the objective function (f)
702 : !> \param last_f the last value of f
703 : !> \param projected_gradient the infinity norm of the projected gradient
704 : !> \par History
705 : !> none
706 : !> \author Fawzi Mohamed
707 : !> @version 2.2002
708 : !> \note
709 : !> private routine
710 : ! **************************************************************************************************
711 2982 : SUBROUTINE cp_opt_gopt_bcast_res(optimizer, n_iter, f, last_f, &
712 : projected_gradient)
713 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
714 : INTEGER, INTENT(out), OPTIONAL :: n_iter
715 : REAL(kind=dp), INTENT(inout), OPTIONAL :: f, last_f, projected_gradient
716 :
717 : REAL(kind=dp), DIMENSION(4) :: results
718 :
719 2982 : IF (optimizer%master == optimizer%para_env%mepos) THEN
720 : results = [REAL(optimizer%isave(30), kind=dp), &
721 7455 : optimizer%f, optimizer%dsave(2), optimizer%dsave(13)]
722 : END IF
723 2982 : CALL optimizer%para_env%bcast(results, optimizer%master)
724 2982 : IF (PRESENT(n_iter)) n_iter = NINT(results(1))
725 2982 : IF (PRESENT(f)) f = results(2)
726 2982 : IF (PRESENT(last_f)) last_f = results(3)
727 2982 : IF (PRESENT(projected_gradient)) projected_gradient = results(4)
728 :
729 2982 : END SUBROUTINE cp_opt_gopt_bcast_res
730 :
731 : ! **************************************************************************************************
732 : !> \brief goes to the next optimal point (after an optimizer iteration)
733 : !> returns true if converged
734 : !> \param optimizer the optimizer that goes to the next point
735 : !> \param n_iter ...
736 : !> \param f ...
737 : !> \param last_f ...
738 : !> \param projected_gradient ...
739 : !> \param converged ...
740 : !> \param geo_section ...
741 : !> \param force_env ...
742 : !> \param gopt_param ...
743 : !> \param spgr ...
744 : !> \return ...
745 : !> \par History
746 : !> 01.2020 modified [pcazade]
747 : !> \author Fawzi Mohamed
748 : !> @version 2.2002
749 : !> \note
750 : !> if you deactivate convergence control it returns never false
751 : ! **************************************************************************************************
752 2982 : FUNCTION cp_opt_gopt_next(optimizer, n_iter, f, last_f, &
753 : projected_gradient, converged, geo_section, force_env, &
754 : gopt_param, spgr) RESULT(res)
755 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
756 : INTEGER, INTENT(out), OPTIONAL :: n_iter
757 : REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
758 : LOGICAL, INTENT(out) :: converged
759 : TYPE(section_vals_type), POINTER :: geo_section
760 : TYPE(force_env_type), POINTER :: force_env
761 : TYPE(gopt_param_type), POINTER :: gopt_param
762 : TYPE(spgr_type), OPTIONAL, POINTER :: spgr
763 : LOGICAL :: res
764 :
765 : ! passes spgr structure if present
766 : CALL cp_opt_gopt_step(optimizer, n_iter=n_iter, f=f, &
767 : last_f=last_f, projected_gradient=projected_gradient, &
768 : converged=converged, geo_section=geo_section, &
769 2982 : force_env=force_env, gopt_param=gopt_param, spgr=spgr)
770 2982 : res = (optimizer%status < 4) .AND. .NOT. converged
771 :
772 2982 : END FUNCTION cp_opt_gopt_next
773 :
774 : ! **************************************************************************************************
775 : !> \brief stops the optimization
776 : !> \param optimizer ...
777 : !> \par History
778 : !> none
779 : !> \author Fawzi Mohamed
780 : !> @version 2.2002
781 : ! **************************************************************************************************
782 0 : SUBROUTINE cp_opt_gopt_stop(optimizer)
783 : TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
784 :
785 0 : optimizer%task = 'STOPPED on user request'
786 0 : optimizer%status = 4 ! normal exit
787 0 : IF (optimizer%master == optimizer%para_env%mepos) THEN
788 0 : CALL cp_opt_gopt_setulb(optimizer)
789 : END IF
790 :
791 0 : END SUBROUTINE cp_opt_gopt_stop
792 :
793 0 : END MODULE cp_lbfgs_optimizer_gopt
|
