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