(git:3add494)
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-2024 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 
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 ! **************************************************************************************************
62 INTERFACE
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 
82 END 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 
177 CONTAINS
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%wanted_relative_f_delta = wanted_relative_f_delta
263  optimizer%wanted_projected_gradient = wanted_projected_gradient
264  optimizer%kind_of_bound = 0
265  IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound
266  IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound
267  IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound
268  IF (PRESENT(trust_radius)) optimizer%trust_radius = trust_radius
269 
270  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
271  optimizer%lower_bound, optimizer%upper_bound, &
272  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
273  optimizer%wanted_relative_f_delta, &
274  optimizer%wanted_projected_gradient, optimizer%work_array, &
275  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
276  optimizer%csave, optimizer%lsave, optimizer%isave, &
277  optimizer%dsave, optimizer%trust_radius)
278  ELSE
279  NULLIFY ( &
280  optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, &
281  optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, &
282  optimizer%dsave, optimizer%work_array)
283  ALLOCATE (optimizer%x(n))
284  optimizer%x(:) = 0.0_dp
285  ALLOCATE (optimizer%gradient(n))
286  optimizer%gradient(:) = 0.0_dp
287  END IF
288  CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
289  optimizer%status = 0
290 
291  CALL timestop(handle)
292 
293  END SUBROUTINE cp_opt_gopt_create
294 
295 ! **************************************************************************************************
296 !> \brief releases the optimizer (see doc/ReferenceCounting.html)
297 !> \param optimizer the object that should be released
298 !> \par History
299 !> 02.2002 created [fawzi]
300 !> 09.2003 dealloc_ref->release [fawzi]
301 !> \author Fawzi Mohamed
302 ! **************************************************************************************************
303  SUBROUTINE cp_opt_gopt_release(optimizer)
304  TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
305 
306  CHARACTER(len=*), PARAMETER :: routinen = 'cp_opt_gopt_release'
307 
308  INTEGER :: handle
309 
310  CALL timeset(routinen, handle)
311 
312  IF (ASSOCIATED(optimizer%kind_of_bound)) THEN
313  DEALLOCATE (optimizer%kind_of_bound)
314  END IF
315  IF (ASSOCIATED(optimizer%i_work_array)) THEN
316  DEALLOCATE (optimizer%i_work_array)
317  END IF
318  IF (ASSOCIATED(optimizer%isave)) THEN
319  DEALLOCATE (optimizer%isave)
320  END IF
321  IF (ASSOCIATED(optimizer%x)) THEN
322  DEALLOCATE (optimizer%x)
323  END IF
324  IF (ASSOCIATED(optimizer%lower_bound)) THEN
325  DEALLOCATE (optimizer%lower_bound)
326  END IF
327  IF (ASSOCIATED(optimizer%upper_bound)) THEN
328  DEALLOCATE (optimizer%upper_bound)
329  END IF
330  IF (ASSOCIATED(optimizer%gradient)) THEN
331  DEALLOCATE (optimizer%gradient)
332  END IF
333  IF (ASSOCIATED(optimizer%dsave)) THEN
334  DEALLOCATE (optimizer%dsave)
335  END IF
336  IF (ASSOCIATED(optimizer%work_array)) THEN
337  DEALLOCATE (optimizer%work_array)
338  END IF
339  CALL mp_para_env_release(optimizer%para_env)
340  CALL gopt_f_release(optimizer%obj_funct)
341 
342  CALL timestop(handle)
343  END SUBROUTINE cp_opt_gopt_release
344 
345 ! **************************************************************************************************
346 !> \brief takes different valuse from the optimizer
347 !> \param optimizer ...
348 !> \param para_env ...
349 !> \param obj_funct ...
350 !> \param m ...
351 !> \param print_every ...
352 !> \param wanted_relative_f_delta ...
353 !> \param wanted_projected_gradient ...
354 !> \param x ...
355 !> \param lower_bound ...
356 !> \param upper_bound ...
357 !> \param kind_of_bound ...
358 !> \param master ...
359 !> \param actual_projected_gradient ...
360 !> \param n_var ...
361 !> \param n_iter ...
362 !> \param status ...
363 !> \param max_f_per_iter ...
364 !> \param at_end ...
365 !> \param is_master ...
366 !> \param last_f ...
367 !> \param f ...
368 !> \par History
369 !> none
370 !> \author Fawzi Mohamed
371 !> @version 2.2002
372 ! **************************************************************************************************
373  SUBROUTINE cp_opt_gopt_get(optimizer, para_env, &
374  obj_funct, m, print_every, &
375  wanted_relative_f_delta, wanted_projected_gradient, &
376  x, lower_bound, upper_bound, kind_of_bound, master, &
377  actual_projected_gradient, &
378  n_var, n_iter, status, max_f_per_iter, at_end, &
379  is_master, last_f, f)
380  TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
381  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
382  TYPE(gopt_f_type), OPTIONAL, POINTER :: obj_funct
383  INTEGER, INTENT(out), OPTIONAL :: m, print_every
384  REAL(kind=dp), INTENT(out), OPTIONAL :: wanted_relative_f_delta, &
385  wanted_projected_gradient
386  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: x, lower_bound, upper_bound
387  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: kind_of_bound
388  INTEGER, INTENT(out), OPTIONAL :: master
389  REAL(kind=dp), INTENT(out), OPTIONAL :: actual_projected_gradient
390  INTEGER, INTENT(out), OPTIONAL :: n_var, n_iter, status, max_f_per_iter
391  LOGICAL, INTENT(out), OPTIONAL :: at_end, is_master
392  REAL(kind=dp), INTENT(out), OPTIONAL :: last_f, f
393 
394  IF (PRESENT(is_master)) is_master = optimizer%master == optimizer%para_env%mepos
395  IF (PRESENT(master)) master = optimizer%master
396  IF (PRESENT(status)) status = optimizer%status
397  IF (PRESENT(para_env)) para_env => optimizer%para_env
398  IF (PRESENT(obj_funct)) obj_funct = optimizer%obj_funct
399  IF (PRESENT(m)) m = optimizer%m
400  IF (PRESENT(max_f_per_iter)) max_f_per_iter = optimizer%max_f_per_iter
401  IF (PRESENT(wanted_projected_gradient)) &
402  wanted_projected_gradient = optimizer%wanted_projected_gradient
403  IF (PRESENT(wanted_relative_f_delta)) &
404  wanted_relative_f_delta = optimizer%wanted_relative_f_delta
405  IF (PRESENT(print_every)) print_every = optimizer%print_every
406  IF (PRESENT(x)) x => optimizer%x
407  IF (PRESENT(n_var)) n_var = SIZE(x)
408  IF (PRESENT(lower_bound)) lower_bound => optimizer%lower_bound
409  IF (PRESENT(upper_bound)) upper_bound => optimizer%upper_bound
410  IF (PRESENT(kind_of_bound)) kind_of_bound => optimizer%kind_of_bound
411  IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
412  IF (PRESENT(last_f)) last_f = optimizer%last_f
413  IF (PRESENT(f)) f = optimizer%f
414  IF (PRESENT(at_end)) at_end = optimizer%status > 3
415  IF (PRESENT(actual_projected_gradient)) &
416  actual_projected_gradient = optimizer%projected_gradient
417  IF (optimizer%master == optimizer%para_env%mepos) THEN
418  IF (optimizer%isave(30) > 1 .AND. (optimizer%task(1:5) == "NEW_X" .OR. &
419  optimizer%task(1:4) == "STOP" .AND. optimizer%task(7:9) == "CPU")) THEN
420  ! nr iterations >1 .and. dsave contains the wanted data
421  IF (PRESENT(last_f)) last_f = optimizer%dsave(2)
422  IF (PRESENT(actual_projected_gradient)) &
423  actual_projected_gradient = optimizer%dsave(13)
424  ELSE
425  cpassert(.NOT. PRESENT(last_f))
426  cpassert(.NOT. PRESENT(actual_projected_gradient))
427  END IF
428  ELSE
429  IF (PRESENT(lower_bound) .OR. PRESENT(upper_bound) .OR. PRESENT(kind_of_bound)) &
430  cpwarn("asked undefined types")
431  END IF
432 
433  END SUBROUTINE cp_opt_gopt_get
434 
435 ! **************************************************************************************************
436 !> \brief does one optimization step
437 !> \param optimizer ...
438 !> \param n_iter ...
439 !> \param f ...
440 !> \param last_f ...
441 !> \param projected_gradient ...
442 !> \param converged ...
443 !> \param geo_section ...
444 !> \param force_env ...
445 !> \param gopt_param ...
446 !> \param spgr ...
447 !> \par History
448 !> 01.2020 modified [pcazade]
449 !> \author Fawzi Mohamed
450 !> @version 2.2002
451 !> \note
452 !> use directly mainlb in place of setulb ??
453 ! **************************************************************************************************
454  SUBROUTINE cp_opt_gopt_step(optimizer, n_iter, f, last_f, &
455  projected_gradient, converged, geo_section, force_env, &
456  gopt_param, spgr)
457  TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
458  INTEGER, INTENT(out), OPTIONAL :: n_iter
459  REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
460  LOGICAL, INTENT(out), OPTIONAL :: converged
461  TYPE(section_vals_type), POINTER :: geo_section
462  TYPE(force_env_type), POINTER :: force_env
463  TYPE(gopt_param_type), POINTER :: gopt_param
464  TYPE(spgr_type), OPTIONAL, POINTER :: spgr
465 
466  CHARACTER(len=*), PARAMETER :: routinen = 'cp_opt_gopt_step'
467 
468  CHARACTER(LEN=5) :: wildcard
469  INTEGER :: dataunit, handle, its
470  LOGICAL :: conv, is_master, justentred, &
471  keep_space_group
472  REAL(kind=dp) :: t_diff, t_now, t_old
473  REAL(kind=dp), DIMENSION(:), POINTER :: xold
474  TYPE(cp_logger_type), POINTER :: logger
475  TYPE(cp_subsys_type), POINTER :: subsys
476 
477  NULLIFY (logger, xold)
478  logger => cp_get_default_logger()
479  CALL timeset(routinen, handle)
480  justentred = .true.
481  is_master = optimizer%master == optimizer%para_env%mepos
482  IF (PRESENT(converged)) converged = optimizer%status == 4
483  ALLOCATE (xold(SIZE(optimizer%x)))
484 
485  ! collecting subsys
486  CALL force_env_get(force_env, subsys=subsys)
487 
488  keep_space_group = .false.
489  IF (PRESENT(spgr)) THEN
490  IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
491  END IF
492 
493  ! applies rotation matrices to coordinates
494  IF (keep_space_group) THEN
495  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
496  END IF
497 
498  xold = optimizer%x
499  t_old = m_walltime()
500 
501  IF (optimizer%status >= 4) THEN
502  cpwarn("status>=4, trying to restart")
503  optimizer%status = 0
504  IF (is_master) THEN
505  optimizer%task = 'START'
506  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
507  optimizer%lower_bound, optimizer%upper_bound, &
508  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
509  optimizer%wanted_relative_f_delta, &
510  optimizer%wanted_projected_gradient, optimizer%work_array, &
511  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
512  optimizer%csave, optimizer%lsave, optimizer%isave, &
513  optimizer%dsave, optimizer%trust_radius, spgr=spgr)
514  END IF
515  END IF
516 
517  DO
518  ifmaster: IF (is_master) THEN
519  IF (optimizer%task(1:7) == 'RESTART') THEN
520  ! restart the optimizer
521  optimizer%status = 0
522  optimizer%task = 'START'
523  ! applies rotation matrices to coordinates and forces
524  IF (keep_space_group) THEN
525  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
526  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
527  END IF
528  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
529  optimizer%lower_bound, optimizer%upper_bound, &
530  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
531  optimizer%wanted_relative_f_delta, &
532  optimizer%wanted_projected_gradient, optimizer%work_array, &
533  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
534  optimizer%csave, optimizer%lsave, optimizer%isave, &
535  optimizer%dsave, optimizer%trust_radius, spgr=spgr)
536  IF (keep_space_group) THEN
537  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
538  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
539  END IF
540  END IF
541  IF (optimizer%task(1:2) == 'FG') THEN
542  IF (optimizer%isave(36) > optimizer%max_f_per_iter) THEN
543  optimizer%task = 'STOP: CPU, hit max f eval in iter'
544  optimizer%status = 5 ! anormal exit
545  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
546  optimizer%lower_bound, optimizer%upper_bound, &
547  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
548  optimizer%wanted_relative_f_delta, &
549  optimizer%wanted_projected_gradient, optimizer%work_array, &
550  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
551  optimizer%csave, optimizer%lsave, optimizer%isave, &
552  optimizer%dsave, optimizer%trust_radius, spgr=spgr)
553  ELSE
554  optimizer%status = 1
555  END IF
556  ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN
557  IF (justentred) THEN
558  optimizer%status = 2
559  ! applies rotation matrices to coordinates and forces
560  IF (keep_space_group) THEN
561  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
562  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
563  END IF
564  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
565  optimizer%lower_bound, optimizer%upper_bound, &
566  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
567  optimizer%wanted_relative_f_delta, &
568  optimizer%wanted_projected_gradient, optimizer%work_array, &
569  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
570  optimizer%csave, optimizer%lsave, optimizer%isave, &
571  optimizer%dsave, optimizer%trust_radius, spgr=spgr)
572  IF (keep_space_group) THEN
573  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
574  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
575  END IF
576  ELSE
577  ! applies rotation matrices to coordinates and forces
578  IF (keep_space_group) THEN
579  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
580  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
581  END IF
582  optimizer%status = 3
583  END IF
584  ELSE IF (optimizer%task(1:4) == 'CONV') THEN
585  optimizer%status = 4
586  ELSE IF (optimizer%task(1:4) == 'STOP') THEN
587  optimizer%status = 5
588  cpwarn("task became stop in an unknown way")
589  ELSE IF (optimizer%task(1:5) == 'ERROR') THEN
590  optimizer%status = 5
591  ELSE
592  cpwarn("unknown task '"//optimizer%task//"'")
593  END IF
594  END IF ifmaster
595  CALL optimizer%para_env%bcast(optimizer%status, optimizer%master)
596  ! Dump info
597  IF (optimizer%status == 3) THEN
598  its = 0
599  IF (is_master) THEN
600  ! Iteration level is taken into account in the optimizer external loop
601  its = optimizer%isave(30)
602  END IF
603  END IF
604  !
605  SELECT CASE (optimizer%status)
606  CASE (1)
607  !op=1 evaluate f and g
608  CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, &
609  f=optimizer%f, &
610  gradient=optimizer%gradient, &
611  final_evaluation=.false., &
612  master=optimizer%master, para_env=optimizer%para_env)
613  ! do not use keywords?
614  IF (is_master) THEN
615  ! applies rotation matrices to coordinates and forces
616  IF (keep_space_group) THEN
617  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
618  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
619  END IF
620  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
621  optimizer%lower_bound, optimizer%upper_bound, &
622  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
623  optimizer%wanted_relative_f_delta, &
624  optimizer%wanted_projected_gradient, optimizer%work_array, &
625  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
626  optimizer%csave, optimizer%lsave, optimizer%isave, &
627  optimizer%dsave, optimizer%trust_radius, spgr=spgr)
628  IF (keep_space_group) THEN
629  CALL spgr_apply_rotations_coord(spgr, optimizer%x)
630  CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
631  END IF
632  END IF
633  CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
634  CASE (2)
635  !op=2 begin new iter
636  CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
637  t_old = m_walltime()
638  CASE (3)
639  !op=3 ended iter
640  wildcard = "LBFGS"
641  dataunit = cp_print_key_unit_nr(logger, geo_section, &
642  "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
643  IF (is_master) its = optimizer%isave(30)
644  CALL optimizer%para_env%bcast(its, optimizer%master)
645 
646  ! Some IO and Convergence check
647  t_now = m_walltime()
648  t_diff = t_now - t_old
649  t_old = t_now
650  CALL gopt_f_io(optimizer%obj_funct, force_env, force_env%root_section, &
651  its, optimizer%f, dataunit, optimizer%eold, optimizer%emin, wildcard, gopt_param, &
652  SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff)
653  CALL optimizer%para_env%bcast(conv, optimizer%master)
654  CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
655  "PRINT%PROGRAM_RUN_INFO")
656  optimizer%eold = optimizer%f
657  optimizer%emin = min(optimizer%emin, optimizer%eold)
658  xold = optimizer%x
659  IF (PRESENT(converged)) converged = conv
660  EXIT
661  CASE (4)
662  !op=4 (convergence - normal exit)
663  ! Specific L-BFGS convergence criteria.. overrides the convergence criteria on
664  ! stepsize and gradients
665  dataunit = cp_print_key_unit_nr(logger, geo_section, &
666  "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
667  IF (dataunit > 0) THEN
668  WRITE (dataunit, '(T2,A)') ""
669  WRITE (dataunit, '(T2,A)') "***********************************************"
670  WRITE (dataunit, '(T2,A)') "* Specific L-BFGS convergence criteria "
671  WRITE (dataunit, '(T2,A)') "* WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR "
672  WRITE (dataunit, '(T2,A)') "* satisfied .... run CONVERGED! "
673  WRITE (dataunit, '(T2,A)') "***********************************************"
674  WRITE (dataunit, '(T2,A)') ""
675  END IF
676  CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
677  "PRINT%PROGRAM_RUN_INFO")
678  IF (PRESENT(converged)) converged = .true.
679  EXIT
680  CASE (5)
681  ! op=5 abnormal exit ()
682  CALL optimizer%para_env%bcast(optimizer%task, optimizer%master)
683  CASE (6)
684  ! deallocated
685  cpabort("step on a deallocated opt structure ")
686  CASE default
687  CALL cp_abort(__location__, &
688  "unknown status "//cp_to_string(optimizer%status))
689  optimizer%status = 5
690  EXIT
691  END SELECT
692  IF (optimizer%status == 1 .AND. justentred) THEN
693  optimizer%eold = optimizer%f
694  optimizer%emin = optimizer%eold
695  END IF
696  justentred = .false.
697  END DO
698 
699  CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
700  CALL cp_opt_gopt_bcast_res(optimizer, &
701  n_iter=optimizer%n_iter, &
702  f=optimizer%f, last_f=optimizer%last_f, &
703  projected_gradient=optimizer%projected_gradient)
704 
705  DEALLOCATE (xold)
706  IF (PRESENT(f)) f = optimizer%f
707  IF (PRESENT(last_f)) last_f = optimizer%last_f
708  IF (PRESENT(projected_gradient)) &
709  projected_gradient = optimizer%projected_gradient
710  IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
711  CALL timestop(handle)
712 
713  END SUBROUTINE cp_opt_gopt_step
714 
715 ! **************************************************************************************************
716 !> \brief returns the results (and broadcasts them)
717 !> \param optimizer the optimizer object the info is taken from
718 !> \param n_iter the number of iterations
719 !> \param f the actual value of the objective function (f)
720 !> \param last_f the last value of f
721 !> \param projected_gradient the infinity norm of the projected gradient
722 !> \par History
723 !> none
724 !> \author Fawzi Mohamed
725 !> @version 2.2002
726 !> \note
727 !> private routine
728 ! **************************************************************************************************
729  SUBROUTINE cp_opt_gopt_bcast_res(optimizer, n_iter, f, last_f, &
730  projected_gradient)
731  TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN) :: optimizer
732  INTEGER, INTENT(out), OPTIONAL :: n_iter
733  REAL(kind=dp), INTENT(inout), OPTIONAL :: f, last_f, projected_gradient
734 
735  REAL(kind=dp), DIMENSION(4) :: results
736 
737  IF (optimizer%master == optimizer%para_env%mepos) THEN
738  results = (/real(optimizer%isave(30), kind=dp), &
739  optimizer%f, optimizer%dsave(2), optimizer%dsave(13)/)
740  END IF
741  CALL optimizer%para_env%bcast(results, optimizer%master)
742  IF (PRESENT(n_iter)) n_iter = nint(results(1))
743  IF (PRESENT(f)) f = results(2)
744  IF (PRESENT(last_f)) last_f = results(3)
745  IF (PRESENT(projected_gradient)) &
746  projected_gradient = results(4)
747 
748  END SUBROUTINE cp_opt_gopt_bcast_res
749 
750 ! **************************************************************************************************
751 !> \brief goes to the next optimal point (after an optimizer iteration)
752 !> returns true if converged
753 !> \param optimizer the optimizer that goes to the next point
754 !> \param n_iter ...
755 !> \param f ...
756 !> \param last_f ...
757 !> \param projected_gradient ...
758 !> \param converged ...
759 !> \param geo_section ...
760 !> \param force_env ...
761 !> \param gopt_param ...
762 !> \param spgr ...
763 !> \return ...
764 !> \par History
765 !> 01.2020 modified [pcazade]
766 !> \author Fawzi Mohamed
767 !> @version 2.2002
768 !> \note
769 !> if you deactivate convergence control it returns never false
770 ! **************************************************************************************************
771  FUNCTION cp_opt_gopt_next(optimizer, n_iter, f, last_f, &
772  projected_gradient, converged, geo_section, force_env, &
773  gopt_param, spgr) RESULT(res)
774  TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
775  INTEGER, INTENT(out), OPTIONAL :: n_iter
776  REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient
777  LOGICAL, INTENT(out) :: converged
778  TYPE(section_vals_type), POINTER :: geo_section
779  TYPE(force_env_type), POINTER :: force_env
780  TYPE(gopt_param_type), POINTER :: gopt_param
781  TYPE(spgr_type), OPTIONAL, POINTER :: spgr
782  LOGICAL :: res
783 
784  ! passes spgr structure if present
785  CALL cp_opt_gopt_step(optimizer, n_iter=n_iter, f=f, &
786  last_f=last_f, projected_gradient=projected_gradient, &
787  converged=converged, geo_section=geo_section, &
788  force_env=force_env, gopt_param=gopt_param, spgr=spgr)
789  res = (optimizer%status < 40) .AND. .NOT. converged
790 
791  END FUNCTION cp_opt_gopt_next
792 
793 ! **************************************************************************************************
794 !> \brief stops the optimization
795 !> \param optimizer ...
796 !> \par History
797 !> none
798 !> \author Fawzi Mohamed
799 !> @version 2.2002
800 ! **************************************************************************************************
801  SUBROUTINE cp_opt_gopt_stop(optimizer)
802  TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT) :: optimizer
803 
804  optimizer%task = 'STOPPED on user request'
805  optimizer%status = 4 ! normal exit
806  IF (optimizer%master == optimizer%para_env%mepos) THEN
807  CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
808  optimizer%lower_bound, optimizer%upper_bound, &
809  optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
810  optimizer%wanted_relative_f_delta, &
811  optimizer%wanted_projected_gradient, optimizer%work_array, &
812  optimizer%i_work_array, optimizer%task, optimizer%print_every, &
813  optimizer%csave, optimizer%lsave, optimizer%isave, &
814  optimizer%dsave, optimizer%trust_radius)
815  END IF
816 
817  END SUBROUTINE cp_opt_gopt_stop
818 
819 END 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_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:774
cp_lbfgs_optimizer_gopt::cp_opt_gopt_stop
subroutine, public cp_opt_gopt_stop(optimizer)
stops the optimization
Definition: cp_lbfgs_optimizer_gopt.F:802
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:304
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)
This subroutine partitions the working arrays wa and iwa, and then uses the limited memory BFGS metho...
Definition: cp_lbfgs.F:185
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::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:1013
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)
returns various attributes about the force environment
Definition: force_env_types.F:329
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:276
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:186
gopt_f_types::gopt_f_release
recursive subroutine, public gopt_f_release(gopt_env)
...
Definition: gopt_f_types.F:200
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:123
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:2027
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:618
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:679