(git:21ef868)
Loading...
Searching...
No Matches
cp_lbfgs_optimizer_gopt.F
Go to the documentation of this file.
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! **************************************************************************************************
22MODULE 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
51! **************************************************************************************************
52!> \brief evaluate the potential energy and its gradients using an array
53!> with same dimension as the particle_set
54!> \param gopt_env the geometry optimization environment
55!> \param x the position where the function should be evaluated
56!> \param f the function value
57!> \param gradient the value of its gradient
58!> \par History
59!> none
60!> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
61! **************************************************************************************************
62INTERFACE
63
64 SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
65 final_evaluation, para_env)
66
67 USE message_passing, ONLY: mp_para_env_type
68 USE gopt_f_types, ONLY: gopt_f_type
69 USE kinds, ONLY: dp
70
71 TYPE(gopt_f_type), POINTER :: gopt_env
72 REAL(KIND=dp), DIMENSION(:), POINTER :: x
73 REAL(KIND=dp), INTENT(out), OPTIONAL :: f
74 REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
75 POINTER :: gradient
76 INTEGER, INTENT(IN) :: master
77 LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
78 TYPE(mp_para_env_type), POINTER :: para_env
79
80 END SUBROUTINE cp_eval_at
81
82END INTERFACE
83
84 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
85 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_lbfgs_optimizer_gopt'
86
87 ! types
88 PUBLIC :: cp_lbfgs_opt_gopt_type
89
90 ! core methods
91
92 ! special methos
93
94 ! underlying functions
95 PUBLIC :: cp_opt_gopt_create, cp_opt_gopt_release, &
96 cp_opt_gopt_next, &
97 cp_opt_gopt_stop
98
99! **************************************************************************************************
100!> \brief info for the optimizer (see the description of this module)
101!> \param task the actual task of the optimizer (in the master it is up to
102!> date, in case of error also the minions one get updated.
103!> \param csave internal character string used by the lbfgs optimizer,
104!> meaningful only in the master
105!> \param lsave logical array used by the lbfgs optimizer, updated only
106!> in the master
107!> On exit with task = 'NEW_X', the following information is
108!> available:
109!> lsave(1) = .true. the initial x did not satisfy the bounds;
110!> lsave(2) = .true. the problem contains bounds;
111!> lsave(3) = .true. each variable has upper and lower bounds.
112!> \param ref_count reference count (see doc/ReferenceCounting.html)
113!> \param m the dimension of the subspace used to approximate the second
114!> derivative
115!> \param print_every every how many iterations output should be written.
116!> if 0 only at end, if print_every<0 never
117!> \param master the pid of the master processor
118!> \param max_f_per_iter the maximum number of function evaluations per
119!> iteration
120!> \param status 0: just initialized, 1: f g calculation,
121!> 2: begin new iteration, 3: ended iteration,
122!> 4: normal (converged) exit, 5: abnormal (error) exit,
123!> 6: daellocated
124!> \param n_iter the actual iteration number
125!> \param kind_of_bound an array with 0 (no bound), 1 (lower bound),
126!> 2 (both bounds), 3 (upper bound), to describe the bounds
127!> of every variable
128!> \param i_work_array an integer workarray of dimension 3*n, present only
129!> in the master
130!> \param isave is an INTEGER working array of dimension 44.
131!> On exit with task = 'NEW_X', it contains information that
132!> the user may want to access:
133!> \param isave (30) = the current iteration number;
134!> \param isave (34) = the total number of function and gradient
135!> evaluations;
136!> \param isave (36) = the number of function value or gradient
137!> evaluations in the current iteration;
138!> \param isave (38) = the number of free variables in the current
139!> iteration;
140!> \param isave (39) = the number of active constraints at the current
141!> iteration;
142!> \param f the actual best value of the object function
143!> \param wanted_relative_f_delta the wanted relative error on f
144!> (to be multiplied by epsilon), 0.0 -> no check
145!> \param wanted_projected_gradient the wanted error on the projected
146!> gradient (hessian times the gradient), 0.0 -> no check
147!> \param last_f the value of f in the last iteration
148!> \param projected_gradient the value of the sup norm of the projected
149!> gradient
150!> \param x the actual evaluation point (best one if converged or stopped)
151!> \param lower_bound the lower bounds
152!> \param upper_bound the upper bounds
153!> \param gradient the actual gradient
154!> \param dsave info date for lbfgs (master only)
155!> \param work_array a work array for lbfgs (master only)
156!> \param para_env the parallel environment for this optimizer
157!> \param obj_funct the objective function to be optimized
158!> \par History
159!> none
160!> \author Fawzi Mohamed
161!> @version 2.2002
162! **************************************************************************************************
163 TYPE cp_lbfgs_opt_gopt_type
164 CHARACTER(len=60) :: task = ""
165 CHARACTER(len=60) :: csave = ""
166 LOGICAL :: lsave(4) = .false.
167 INTEGER :: m = 0, print_every = 0, master = 0, max_f_per_iter = 0, status = 0, n_iter = 0
168 INTEGER, DIMENSION(:), POINTER :: kind_of_bound => null(), i_work_array => null(), isave => null()
169 REAL(kind=dp) :: f = 0.0_dp, wanted_relative_f_delta = 0.0_dp, wanted_projected_gradient = 0.0_dp, &
170 last_f = 0.0_dp, projected_gradient = 0.0_dp, eold = 0.0_dp, emin = 0.0_dp, trust_radius = 0.0_dp
171 REAL(kind=dp), DIMENSION(:), POINTER :: x => null(), lower_bound => null(), upper_bound => null(), &
172 gradient => null(), dsave => null(), work_array => null()
173 TYPE(mp_para_env_type), POINTER :: para_env => null()
174 TYPE(gopt_f_type), POINTER :: obj_funct => null()
175 END TYPE cp_lbfgs_opt_gopt_type
176
177CONTAINS
178
179! **************************************************************************************************
180!> \brief initializes the optimizer
181!> \param optimizer ...
182!> \param para_env ...
183!> \param obj_funct ...
184!> \param x0 ...
185!> \param m ...
186!> \param print_every ...
187!> \param wanted_relative_f_delta ...
188!> \param wanted_projected_gradient ...
189!> \param lower_bound ...
190!> \param upper_bound ...
191!> \param kind_of_bound ...
192!> \param master ...
193!> \param max_f_per_iter ...
194!> \param trust_radius ...
195!> \par History
196!> 02.2002 created [fawzi]
197!> 09.2003 refactored (retain/release,para_env) [fawzi]
198!> \author Fawzi Mohamed
199!> \note
200!> redirects the lbfgs output the the default unit
201! **************************************************************************************************
202 SUBROUTINE cp_opt_gopt_create(optimizer, para_env, obj_funct, x0, m, print_every, &
203 wanted_relative_f_delta, wanted_projected_gradient, lower_bound, upper_bound, &
204 kind_of_bound, master, max_f_per_iter, trust_radius)
205 TYPE(cp_lbfgs_opt_gopt_type), INTENT(OUT) :: optimizer
206 TYPE(mp_para_env_type), POINTER :: para_env
207 TYPE(gopt_f_type), POINTER :: obj_funct
208 REAL(kind=dp), DIMENSION(:), INTENT(in) :: x0
209 INTEGER, INTENT(in), OPTIONAL :: m, print_every
210 REAL(kind=dp), INTENT(in), OPTIONAL :: wanted_relative_f_delta, &
211 wanted_projected_gradient
212 REAL(kind=dp), DIMENSION(SIZE(x0)), INTENT(in), &
213 OPTIONAL :: lower_bound, upper_bound
214 INTEGER, DIMENSION(SIZE(x0)), INTENT(in), OPTIONAL :: kind_of_bound
215 INTEGER, INTENT(in), OPTIONAL :: master, max_f_per_iter
216 REAL(kind=dp), INTENT(in), OPTIONAL :: trust_radius
217
218 CHARACTER(len=*), PARAMETER :: routinen = 'cp_opt_gopt_create'
219
220 INTEGER :: handle, lenwa, n
221
222 CALL timeset(routinen, handle)
223
224 NULLIFY (optimizer%kind_of_bound, &
225 optimizer%i_work_array, &
226 optimizer%isave, &
227 optimizer%x, &
228 optimizer%lower_bound, &
229 optimizer%upper_bound, &
230 optimizer%gradient, &
231 optimizer%dsave, &
232 optimizer%work_array, &
233 optimizer%para_env, &
234 optimizer%obj_funct)
235 n = SIZE(x0)
236 optimizer%m = 4
237 IF (PRESENT(m)) optimizer%m = m
238 optimizer%master = para_env%source
239 optimizer%para_env => para_env
240 CALL para_env%retain()
241 optimizer%obj_funct => obj_funct
242 CALL gopt_f_retain(obj_funct)
243 optimizer%max_f_per_iter = 20
244 IF (PRESENT(max_f_per_iter)) optimizer%max_f_per_iter = max_f_per_iter
245 optimizer%print_every = -1
246 optimizer%n_iter = 0
247 optimizer%f = -1.0_dp
248 optimizer%last_f = -1.0_dp
249 optimizer%projected_gradient = -1.0_dp
250 IF (PRESENT(print_every)) optimizer%print_every = print_every
251 IF (PRESENT(master)) optimizer%master = master
252 IF (optimizer%master == optimizer%para_env%mepos) THEN
253 !MK This has to be adapted for a new L-BFGS version possibly
254 lenwa = 2*optimizer%m*n + 5*n + 11*optimizer%m*optimizer%m + 8*optimizer%m
255 ALLOCATE (optimizer%kind_of_bound(n), optimizer%i_work_array(3*n), &
256 optimizer%isave(44))
257 ALLOCATE (optimizer%x(n), optimizer%lower_bound(n), &
258 optimizer%upper_bound(n), optimizer%gradient(n), &
259 optimizer%dsave(29), optimizer%work_array(lenwa))
260 optimizer%x = x0
261 optimizer%task = 'START'
262 optimizer%i_work_array = 0
263 optimizer%isave = 0
264 optimizer%lower_bound = 0.0_dp
265 optimizer%upper_bound = 0.0_dp
266 optimizer%gradient = 0.0_dp
267 optimizer%dsave = 0.0_dp
268 optimizer%work_array = 0.0_dp
269 IF (PRESENT(wanted_relative_f_delta)) THEN
270 optimizer%wanted_relative_f_delta = wanted_relative_f_delta
271 END IF
272 IF (PRESENT(wanted_projected_gradient)) THEN
273 optimizer%wanted_projected_gradient = wanted_projected_gradient
274 END IF
275 optimizer%kind_of_bound = 0
276 IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound
277 IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound
278 IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound
279 IF (PRESENT(trust_radius)) optimizer%trust_radius = trust_radius
280
281 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
282 optimizer%lower_bound, optimizer%upper_bound, &
283 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
284 optimizer%wanted_relative_f_delta, &
285 optimizer%wanted_projected_gradient, optimizer%work_array, &
286 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
287 optimizer%csave, optimizer%lsave, optimizer%isave, &
288 optimizer%dsave, optimizer%trust_radius)
289 ELSE
290 NULLIFY ( &
291 optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, &
292 optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, &
293 optimizer%dsave, optimizer%work_array)
294 ALLOCATE (optimizer%x(n))
295 optimizer%x(:) = 0.0_dp
296 ALLOCATE (optimizer%gradient(n))
297 optimizer%gradient(:) = 0.0_dp
298 END IF
299 CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
300 optimizer%status = 0
301
302 CALL timestop(handle)
303
304 END SUBROUTINE cp_opt_gopt_create
305
306! **************************************************************************************************
307!> \brief releases the optimizer (see doc/ReferenceCounting.html)
308!> \param optimizer the object that should be released
309!> \par History
310!> 02.2002 created [fawzi]
311!> 09.2003 dealloc_ref->release [fawzi]
312!> \author Fawzi Mohamed
313! **************************************************************************************************
314 SUBROUTINE cp_opt_gopt_release(optimizer)
315 TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
316
317 CHARACTER(len=*), PARAMETER :: routinen = 'cp_opt_gopt_release'
318
319 INTEGER :: handle
320
321 CALL timeset(routinen, handle)
322
323 IF (ASSOCIATED(optimizer%kind_of_bound)) THEN
324 DEALLOCATE (optimizer%kind_of_bound)
325 END IF
326 IF (ASSOCIATED(optimizer%i_work_array)) THEN
327 DEALLOCATE (optimizer%i_work_array)
328 END IF
329 IF (ASSOCIATED(optimizer%isave)) THEN
330 DEALLOCATE (optimizer%isave)
331 END IF
332 IF (ASSOCIATED(optimizer%x)) THEN
333 DEALLOCATE (optimizer%x)
334 END IF
335 IF (ASSOCIATED(optimizer%lower_bound)) THEN
336 DEALLOCATE (optimizer%lower_bound)
337 END IF
338 IF (ASSOCIATED(optimizer%upper_bound)) THEN
339 DEALLOCATE (optimizer%upper_bound)
340 END IF
341 IF (ASSOCIATED(optimizer%gradient)) THEN
342 DEALLOCATE (optimizer%gradient)
343 END IF
344 IF (ASSOCIATED(optimizer%dsave)) THEN
345 DEALLOCATE (optimizer%dsave)
346 END IF
347 IF (ASSOCIATED(optimizer%work_array)) THEN
348 DEALLOCATE (optimizer%work_array)
349 END IF
350 CALL mp_para_env_release(optimizer%para_env)
351 CALL gopt_f_release(optimizer%obj_funct)
352
353 CALL timestop(handle)
354 END SUBROUTINE cp_opt_gopt_release
355
356! **************************************************************************************************
357!> \brief takes different valuse from the optimizer
358!> \param optimizer ...
359!> \param para_env ...
360!> \param obj_funct ...
361!> \param m ...
362!> \param print_every ...
363!> \param wanted_relative_f_delta ...
364!> \param wanted_projected_gradient ...
365!> \param x ...
366!> \param lower_bound ...
367!> \param upper_bound ...
368!> \param kind_of_bound ...
369!> \param master ...
370!> \param actual_projected_gradient ...
371!> \param n_var ...
372!> \param n_iter ...
373!> \param status ...
374!> \param max_f_per_iter ...
375!> \param at_end ...
376!> \param is_master ...
377!> \param last_f ...
378!> \param f ...
379!> \par History
380!> none
381!> \author Fawzi Mohamed
382!> @version 2.2002
383! **************************************************************************************************
384 SUBROUTINE cp_opt_gopt_get(optimizer, para_env, &
385 obj_funct, m, print_every, &
386 wanted_relative_f_delta, wanted_projected_gradient, &
387 x, lower_bound, upper_bound, kind_of_bound, master, &
388 actual_projected_gradient, &
389 n_var, n_iter, status, max_f_per_iter, at_end, &
390 is_master, last_f, f)
391 TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
392 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
393 TYPE(gopt_f_type), OPTIONAL, POINTER :: obj_funct
394 INTEGER, INTENT(out), OPTIONAL :: m, print_every
395 REAL(kind=dp), INTENT(out), OPTIONAL :: wanted_relative_f_delta, &
396 wanted_projected_gradient
397 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: x, lower_bound, upper_bound
398 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: kind_of_bound
399 INTEGER, INTENT(out), OPTIONAL :: master
400 REAL(kind=dp), INTENT(out), OPTIONAL :: actual_projected_gradient
401 INTEGER, INTENT(out), OPTIONAL :: n_var, n_iter, status, max_f_per_iter
402 LOGICAL, INTENT(out), OPTIONAL :: at_end, is_master
403 REAL(kind=dp), INTENT(out), OPTIONAL :: last_f, f
404
405 IF (PRESENT(is_master)) is_master = optimizer%master == optimizer%para_env%mepos
406 IF (PRESENT(master)) master = optimizer%master
407 IF (PRESENT(status)) status = optimizer%status
408 IF (PRESENT(para_env)) para_env => optimizer%para_env
409 IF (PRESENT(obj_funct)) obj_funct = optimizer%obj_funct
410 IF (PRESENT(m)) m = optimizer%m
411 IF (PRESENT(max_f_per_iter)) max_f_per_iter = optimizer%max_f_per_iter
412 IF (PRESENT(wanted_projected_gradient)) THEN
413 wanted_projected_gradient = optimizer%wanted_projected_gradient
414 END IF
415 IF (PRESENT(wanted_relative_f_delta)) THEN
416 wanted_relative_f_delta = optimizer%wanted_relative_f_delta
417 END IF
418 IF (PRESENT(print_every)) print_every = optimizer%print_every
419 IF (PRESENT(x)) x => optimizer%x
420 IF (PRESENT(n_var)) n_var = SIZE(x)
421 IF (PRESENT(lower_bound)) lower_bound => optimizer%lower_bound
422 IF (PRESENT(upper_bound)) upper_bound => optimizer%upper_bound
423 IF (PRESENT(kind_of_bound)) kind_of_bound => optimizer%kind_of_bound
424 IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
425 IF (PRESENT(last_f)) last_f = optimizer%last_f
426 IF (PRESENT(f)) f = optimizer%f
427 IF (PRESENT(at_end)) at_end = optimizer%status > 3
428 IF (PRESENT(actual_projected_gradient)) THEN
429 actual_projected_gradient = optimizer%projected_gradient
430 END IF
431 IF (optimizer%master == optimizer%para_env%mepos) THEN
432 IF (optimizer%isave(30) > 1 .AND. (optimizer%task(1:5) == "NEW_X" .OR. &
433 optimizer%task(1:4) == "STOP" .AND. optimizer%task(7:9) == "CPU")) THEN
434 ! nr iterations >1 .and. dsave contains the wanted data
435 IF (PRESENT(last_f)) last_f = optimizer%dsave(2)
436 IF (PRESENT(actual_projected_gradient)) THEN
437 actual_projected_gradient = optimizer%dsave(13)
438 END IF
439 ELSE
440 cpassert(.NOT. PRESENT(last_f))
441 cpassert(.NOT. PRESENT(actual_projected_gradient))
442 END IF
443 ELSE IF (PRESENT(lower_bound) .OR. PRESENT(upper_bound) .OR. PRESENT(kind_of_bound)) THEN
444 cpwarn("asked undefined types")
445 END IF
446
447 END SUBROUTINE cp_opt_gopt_get
448
449! **************************************************************************************************
450!> \brief does one optimization step
451!> \param optimizer ...
452!> \param n_iter ...
453!> \param f ...
454!> \param last_f ...
455!> \param projected_gradient ...
456!> \param converged ...
457!> \param geo_section ...
458!> \param force_env ...
459!> \param gopt_param ...
460!> \param spgr ...
461!> \par History
462!> 01.2020 modified [pcazade]
463!> \author Fawzi Mohamed
464!> @version 2.2002
465!> \note
466!> use directly mainlb in place of setulb ??
467! **************************************************************************************************
468 SUBROUTINE cp_opt_gopt_step(optimizer, n_iter, f, last_f, &
469 projected_gradient, converged, geo_section, force_env, &
470 gopt_param, spgr)
471 TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
472 INTEGER, INTENT(out), OPTIONAL :: n_iter
473 REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
474 LOGICAL, INTENT(out), OPTIONAL :: converged
475 TYPE(section_vals_type), POINTER :: geo_section
476 TYPE(force_env_type), POINTER :: force_env
477 TYPE(gopt_param_type), POINTER :: gopt_param
478 TYPE(spgr_type), OPTIONAL, POINTER :: spgr
479
480 CHARACTER(len=*), PARAMETER :: routinen = 'cp_opt_gopt_step'
481
482 CHARACTER(LEN=5) :: wildcard
483 INTEGER :: dataunit, handle, its
484 LOGICAL :: conv, is_master, justentred, &
485 keep_space_group
486 REAL(kind=dp) :: t_diff, t_now, t_old
487 REAL(kind=dp), DIMENSION(:), POINTER :: xold
488 TYPE(cp_logger_type), POINTER :: logger
489 TYPE(cp_subsys_type), POINTER :: subsys
490
491 NULLIFY (logger, xold)
492 logger => cp_get_default_logger()
493 CALL timeset(routinen, handle)
494 justentred = .true.
495 is_master = optimizer%master == optimizer%para_env%mepos
496 IF (PRESENT(converged)) converged = optimizer%status == 4
497 ALLOCATE (xold(SIZE(optimizer%x)))
498
499 ! collecting subsys
500 CALL force_env_get(force_env, subsys=subsys)
501
502 keep_space_group = .false.
503 IF (PRESENT(spgr)) THEN
504 IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
505 END IF
506
507 ! applies rotation matrices to coordinates
508 IF (keep_space_group) THEN
509 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
510 END IF
511
512 xold = optimizer%x
513 t_old = m_walltime()
514
515 IF (optimizer%status >= 4) THEN
516 cpwarn("status>=4, trying to restart")
517 optimizer%status = 0
518 dataunit = cp_print_key_unit_nr(logger, geo_section, &
519 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
520 IF (is_master) THEN
521 optimizer%task = 'START'
522 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
523 optimizer%lower_bound, optimizer%upper_bound, &
524 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
525 optimizer%wanted_relative_f_delta, &
526 optimizer%wanted_projected_gradient, optimizer%work_array, &
527 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
528 optimizer%csave, optimizer%lsave, optimizer%isave, &
529 optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=dataunit)
530 END IF
531 CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
532 "PRINT%PROGRAM_RUN_INFO")
533 END IF
534
535 DO
536 dataunit = cp_print_key_unit_nr(logger, geo_section, &
537 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
538 ifmaster: IF (is_master) THEN
539 IF (optimizer%task(1:7) == 'RESTART') THEN
540 ! restart the optimizer
541 optimizer%status = 0
542 optimizer%task = 'START'
543 ! applies rotation matrices to coordinates and forces
544 IF (keep_space_group) THEN
545 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
546 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
547 END IF
548 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
549 optimizer%lower_bound, optimizer%upper_bound, &
550 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
551 optimizer%wanted_relative_f_delta, &
552 optimizer%wanted_projected_gradient, optimizer%work_array, &
553 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
554 optimizer%csave, optimizer%lsave, optimizer%isave, &
555 optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=dataunit)
556 IF (keep_space_group) THEN
557 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
558 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
559 END IF
560 END IF
561 IF (optimizer%task(1:2) == 'FG') THEN
562 IF (optimizer%isave(36) > optimizer%max_f_per_iter) THEN
563 optimizer%task = 'STOP: CPU, hit max f eval in iter'
564 optimizer%status = 5 ! anormal exit
565 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
566 optimizer%lower_bound, optimizer%upper_bound, &
567 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
568 optimizer%wanted_relative_f_delta, &
569 optimizer%wanted_projected_gradient, optimizer%work_array, &
570 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
571 optimizer%csave, optimizer%lsave, optimizer%isave, &
572 optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=dataunit)
573 ELSE
574 optimizer%status = 1
575 END IF
576 ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN
577 IF (justentred) THEN
578 optimizer%status = 2
579 ! applies rotation matrices to coordinates and forces
580 IF (keep_space_group) THEN
581 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
582 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
583 END IF
584 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
585 optimizer%lower_bound, optimizer%upper_bound, &
586 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
587 optimizer%wanted_relative_f_delta, &
588 optimizer%wanted_projected_gradient, optimizer%work_array, &
589 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
590 optimizer%csave, optimizer%lsave, optimizer%isave, &
591 optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=dataunit)
592 IF (keep_space_group) THEN
593 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
594 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
595 END IF
596 ELSE
597 ! applies rotation matrices to coordinates and forces
598 IF (keep_space_group) THEN
599 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
600 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
601 END IF
602 optimizer%status = 3
603 END IF
604 ELSE IF (optimizer%task(1:4) == 'CONV') THEN
605 optimizer%status = 4
606 ELSE IF (optimizer%task(1:4) == 'STOP') THEN
607 optimizer%status = 5
608 cpwarn("task became stop in an unknown way")
609 ELSE IF (optimizer%task(1:5) == 'ERROR') THEN
610 optimizer%status = 5
611 ELSE
612 cpwarn("unknown task '"//optimizer%task//"'")
613 END IF
614 END IF ifmaster
615 CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
616 "PRINT%PROGRAM_RUN_INFO")
617 CALL optimizer%para_env%bcast(optimizer%status, optimizer%master)
618 ! Dump info
619 IF (optimizer%status == 3) THEN
620 its = 0
621 IF (is_master) THEN
622 ! Iteration level is taken into account in the optimizer external loop
623 its = optimizer%isave(30)
624 END IF
625 END IF
626 !
627 SELECT CASE (optimizer%status)
628 CASE (1)
629 !op=1 evaluate f and g
630 CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, &
631 f=optimizer%f, &
632 gradient=optimizer%gradient, &
633 final_evaluation=.false., &
634 master=optimizer%master, para_env=optimizer%para_env)
635 ! do not use keywords?
636 dataunit = cp_print_key_unit_nr(logger, geo_section, &
637 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
638 IF (is_master) THEN
639 ! applies rotation matrices to coordinates and forces
640 IF (keep_space_group) THEN
641 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
642 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
643 END IF
644 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
645 optimizer%lower_bound, optimizer%upper_bound, &
646 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
647 optimizer%wanted_relative_f_delta, &
648 optimizer%wanted_projected_gradient, optimizer%work_array, &
649 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
650 optimizer%csave, optimizer%lsave, optimizer%isave, &
651 optimizer%dsave, optimizer%trust_radius, spgr=spgr, iwunit=dataunit)
652 IF (keep_space_group) THEN
653 CALL spgr_apply_rotations_coord(spgr, optimizer%x)
654 CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
655 END IF
656 END IF
657 CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
658 "PRINT%PROGRAM_RUN_INFO")
659 CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
660 CASE (2)
661 !op=2 begin new iter
662 CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
663 t_old = m_walltime()
664 CASE (3)
665 !op=3 ended iter
666 wildcard = "LBFGS"
667 dataunit = cp_print_key_unit_nr(logger, geo_section, &
668 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
669 IF (is_master) its = optimizer%isave(30)
670 CALL optimizer%para_env%bcast(its, optimizer%master)
671
672 ! Some IO and Convergence check
673 t_now = m_walltime()
674 t_diff = t_now - t_old
675 t_old = t_now
676 CALL gopt_f_io(optimizer%obj_funct, force_env, force_env%root_section, &
677 its, optimizer%f, dataunit, optimizer%eold, optimizer%emin, wildcard, gopt_param, &
678 SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff)
679 CALL optimizer%para_env%bcast(conv, optimizer%master)
680 CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
681 "PRINT%PROGRAM_RUN_INFO")
682 optimizer%eold = optimizer%f
683 optimizer%emin = min(optimizer%emin, optimizer%eold)
684 xold = optimizer%x
685 IF (PRESENT(converged)) converged = conv
686 EXIT
687 CASE (4)
688 !op=4 (convergence - normal exit)
689 ! Specific L-BFGS convergence criteria.. overrides the convergence criteria on
690 ! stepsize and gradients
691 dataunit = cp_print_key_unit_nr(logger, geo_section, &
692 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
693 IF (dataunit > 0) THEN
694 WRITE (dataunit, '(T2,A)') ""
695 WRITE (dataunit, '(T2,A)') "************************************************"
696 WRITE (dataunit, '(T2,A)') "* Specific L-BFGS convergence criteria *"
697 WRITE (dataunit, '(T2,A)') "* WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR *"
698 WRITE (dataunit, '(T2,A)') "* satisfied .... run CONVERGED! *"
699 WRITE (dataunit, '(T2,A)') "* * * * *"
700 WRITE (dataunit, '(T2,A)') "* General convergence criteria on stepsize and *"
701 WRITE (dataunit, '(T2,A)') "* gradients may or may not have been satisfied *"
702 WRITE (dataunit, '(T2,A)') "* yet; if unsatisfactory, try tightening the *"
703 WRITE (dataunit, '(T2,A)') "* L-BFGS convergence criteria and restart run. *"
704 WRITE (dataunit, '(T2,A)') "************************************************"
705 WRITE (dataunit, '(T2,A)') ""
706 END IF
707 CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
708 "PRINT%PROGRAM_RUN_INFO")
709 IF (PRESENT(converged)) converged = .true.
710 EXIT
711 CASE (5)
712 ! op=5 abnormal exit ()
713 CALL optimizer%para_env%bcast(optimizer%task, optimizer%master)
714 CASE (6)
715 ! deallocated
716 cpabort("step on a deallocated opt structure ")
717 CASE default
718 CALL cp_abort(__location__, &
719 "unknown status "//cp_to_string(optimizer%status))
720 optimizer%status = 5
721 EXIT
722 END SELECT
723 IF (optimizer%status == 1 .AND. justentred) THEN
724 optimizer%eold = optimizer%f
725 optimizer%emin = optimizer%eold
726 END IF
727 justentred = .false.
728 END DO
729
730 CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
731 CALL cp_opt_gopt_bcast_res(optimizer, &
732 n_iter=optimizer%n_iter, &
733 f=optimizer%f, last_f=optimizer%last_f, &
734 projected_gradient=optimizer%projected_gradient)
735
736 DEALLOCATE (xold)
737 IF (PRESENT(f)) f = optimizer%f
738 IF (PRESENT(last_f)) last_f = optimizer%last_f
739 IF (PRESENT(projected_gradient)) projected_gradient = optimizer%projected_gradient
740 IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
741 CALL timestop(handle)
742
743 END SUBROUTINE cp_opt_gopt_step
744
745! **************************************************************************************************
746!> \brief returns the results (and broadcasts them)
747!> \param optimizer the optimizer object the info is taken from
748!> \param n_iter the number of iterations
749!> \param f the actual value of the objective function (f)
750!> \param last_f the last value of f
751!> \param projected_gradient the infinity norm of the projected gradient
752!> \par History
753!> none
754!> \author Fawzi Mohamed
755!> @version 2.2002
756!> \note
757!> private routine
758! **************************************************************************************************
759 SUBROUTINE cp_opt_gopt_bcast_res(optimizer, n_iter, f, last_f, &
760 projected_gradient)
761 TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
762 INTEGER, INTENT(out), OPTIONAL :: n_iter
763 REAL(kind=dp), INTENT(inout), OPTIONAL :: f, last_f, projected_gradient
764
765 REAL(kind=dp), DIMENSION(4) :: results
766
767 IF (optimizer%master == optimizer%para_env%mepos) THEN
768 results = [real(optimizer%isave(30), kind=dp), &
769 optimizer%f, optimizer%dsave(2), optimizer%dsave(13)]
770 END IF
771 CALL optimizer%para_env%bcast(results, optimizer%master)
772 IF (PRESENT(n_iter)) n_iter = nint(results(1))
773 IF (PRESENT(f)) f = results(2)
774 IF (PRESENT(last_f)) last_f = results(3)
775 IF (PRESENT(projected_gradient)) projected_gradient = results(4)
776
777 END SUBROUTINE cp_opt_gopt_bcast_res
778
779! **************************************************************************************************
780!> \brief goes to the next optimal point (after an optimizer iteration)
781!> returns true if converged
782!> \param optimizer the optimizer that goes to the next point
783!> \param n_iter ...
784!> \param f ...
785!> \param last_f ...
786!> \param projected_gradient ...
787!> \param converged ...
788!> \param geo_section ...
789!> \param force_env ...
790!> \param gopt_param ...
791!> \param spgr ...
792!> \return ...
793!> \par History
794!> 01.2020 modified [pcazade]
795!> \author Fawzi Mohamed
796!> @version 2.2002
797!> \note
798!> if you deactivate convergence control it returns never false
799! **************************************************************************************************
800 FUNCTION cp_opt_gopt_next(optimizer, n_iter, f, last_f, &
801 projected_gradient, converged, geo_section, force_env, &
802 gopt_param, spgr) RESULT(res)
803 TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
804 INTEGER, INTENT(out), OPTIONAL :: n_iter
805 REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
806 LOGICAL, INTENT(out) :: converged
807 TYPE(section_vals_type), POINTER :: geo_section
808 TYPE(force_env_type), POINTER :: force_env
809 TYPE(gopt_param_type), POINTER :: gopt_param
810 TYPE(spgr_type), OPTIONAL, POINTER :: spgr
811 LOGICAL :: res
812
813 ! passes spgr structure if present
814 CALL cp_opt_gopt_step(optimizer, n_iter=n_iter, f=f, &
815 last_f=last_f, projected_gradient=projected_gradient, &
816 converged=converged, geo_section=geo_section, &
817 force_env=force_env, gopt_param=gopt_param, spgr=spgr)
818 res = (optimizer%status < 40) .AND. .NOT. converged
819
820 END FUNCTION cp_opt_gopt_next
821
822! **************************************************************************************************
823!> \brief stops the optimization
824!> \param optimizer ...
825!> \par History
826!> none
827!> \author Fawzi Mohamed
828!> @version 2.2002
829! **************************************************************************************************
830 SUBROUTINE cp_opt_gopt_stop(optimizer)
831 TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
832
833 optimizer%task = 'STOPPED on user request'
834 optimizer%status = 4 ! normal exit
835 IF (optimizer%master == optimizer%para_env%mepos) THEN
836 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
837 optimizer%lower_bound, optimizer%upper_bound, &
838 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
839 optimizer%wanted_relative_f_delta, &
840 optimizer%wanted_projected_gradient, optimizer%work_array, &
841 optimizer%i_work_array, optimizer%task, optimizer%print_every, &
842 optimizer%csave, optimizer%lsave, optimizer%isave, &
843 optimizer%dsave, optimizer%trust_radius)
844 END IF
845
846 END SUBROUTINE cp_opt_gopt_stop
847
848END MODULE cp_lbfgs_optimizer_gopt
cp_eval_at
subroutine cp_eval_at(gopt_env, x, f, gradient, master, final_evaluation, para_env)
evaluete the potential energy and its gradients using an array with same dimension as the particle_se...
Definition gopt_f77_methods.F:25
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_lbfgs_optimizer_gopt
routines that optimize a functional using the limited memory bfgs quasi-newton method....
Definition cp_lbfgs_optimizer_gopt.F:22
cp_lbfgs_optimizer_gopt::cp_opt_gopt_create
subroutine, public cp_opt_gopt_create(optimizer, para_env, obj_funct, x0, m, print_every, wanted_relative_f_delta, wanted_projected_gradient, lower_bound, upper_bound, kind_of_bound, master, max_f_per_iter, trust_radius)
initializes the optimizer
Definition cp_lbfgs_optimizer_gopt.F:205
cp_lbfgs_optimizer_gopt::cp_opt_gopt_next
logical function, public cp_opt_gopt_next(optimizer, n_iter, f, last_f, projected_gradient, converged, geo_section, force_env, gopt_param, spgr)
goes to the next optimal point (after an optimizer iteration) returns true if converged
Definition cp_lbfgs_optimizer_gopt.F:803
cp_lbfgs_optimizer_gopt::cp_opt_gopt_stop
subroutine, public cp_opt_gopt_stop(optimizer)
stops the optimization
Definition cp_lbfgs_optimizer_gopt.F:831
cp_lbfgs_optimizer_gopt::cp_opt_gopt_release
subroutine, public cp_opt_gopt_release(optimizer)
releases the optimizer (see doc/ReferenceCounting.html)
Definition cp_lbfgs_optimizer_gopt.F:315
cp_lbfgs
LBFGS-B routine (version 3.0, April 25, 2011)
Definition cp_lbfgs.F:19
cp_lbfgs::setulb
subroutine, public setulb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave, trust_radius, spgr, iwunit)
This subroutine partitions the working arrays wa and iwa, and then uses the limited memory BFGS metho...
Definition cp_lbfgs.F:188
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
gopt_f_methods
contains a functional that calculates the energy and its derivatives for the geometry optimizer
Definition gopt_f_methods.F:14
gopt_f_methods::gopt_f_io
subroutine, public gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, step, rad, used_time)
Handles the Output during an optimization run.
Definition gopt_f_methods.F:315
gopt_f_types
contains a functional that calculates the energy and its derivatives for the geometry optimizer
Definition gopt_f_types.F:15
gopt_f_types::gopt_f_retain
subroutine, public gopt_f_retain(gopt_env)
...
Definition gopt_f_types.F:187
gopt_f_types::gopt_f_release
recursive subroutine, public gopt_f_release(gopt_env)
...
Definition gopt_f_types.F:201
gopt_param_types
contains typo and related routines to handle parameters controlling the GEO_OPT module
Definition gopt_param_types.F:14
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition message_passing.F:2029
space_groups_types
Space Group Symmetry Type Module (version 1.0, Ferbruary 12, 2021)
Definition space_groups_types.F:14
space_groups
Space Group Symmetry Module (version 1.0, January 16, 2020)
Definition space_groups.F:14
space_groups::spgr_apply_rotations_coord
subroutine, public spgr_apply_rotations_coord(spgr, coord)
routine applies the rotation matrices to the coordinates.
Definition space_groups.F:610
space_groups::spgr_apply_rotations_force
subroutine, public spgr_apply_rotations_force(spgr, force)
routine applies the rotation matrices to the forces.
Definition space_groups.F:671
cp_lbfgs_optimizer_gopt::cp_lbfgs_opt_gopt_type
info for the optimizer (see the description of this module)
Definition cp_lbfgs_optimizer_gopt.F:163
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
gopt_f_types::gopt_f_type
calculates the potential energy of a system, and its derivatives
Definition gopt_f_types.F:60
gopt_param_types::gopt_param_type
Definition gopt_param_types.F:55
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
space_groups_types::spgr_type
Definition space_groups_types.F:27