LCOV - code coverage report
Current view: top level - src/motion - cp_lbfgs_optimizer_gopt.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 65.3 % 259 169
Test Date: 2025-07-25 12:55:17 Functions: 62.5 % 8 5

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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         2448 :    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          408 :       CALL timeset(routineN, handle)
     191              : 
     192              :       NULLIFY (optimizer%kind_of_bound, &
     193          408 :                optimizer%i_work_array, &
     194          408 :                optimizer%isave, &
     195          408 :                optimizer%x, &
     196          408 :                optimizer%lower_bound, &
     197          408 :                optimizer%upper_bound, &
     198          408 :                optimizer%gradient, &
     199          408 :                optimizer%dsave, &
     200          408 :                optimizer%work_array, &
     201              :                optimizer%para_env, &
     202          408 :                optimizer%obj_funct)
     203          408 :       n = SIZE(x0)
     204          408 :       optimizer%m = 4
     205          408 :       IF (PRESENT(m)) optimizer%m = m
     206          408 :       optimizer%master = para_env%source
     207          408 :       optimizer%para_env => para_env
     208          408 :       CALL para_env%retain()
     209          408 :       optimizer%obj_funct => obj_funct
     210          408 :       CALL gopt_f_retain(obj_funct)
     211          408 :       optimizer%max_f_per_iter = 20
     212          408 :       IF (PRESENT(max_f_per_iter)) optimizer%max_f_per_iter = max_f_per_iter
     213          408 :       optimizer%print_every = -1
     214          408 :       optimizer%n_iter = 0
     215          408 :       optimizer%f = -1.0_dp
     216          408 :       optimizer%last_f = -1.0_dp
     217          408 :       optimizer%projected_gradient = -1.0_dp
     218          408 :       IF (PRESENT(print_every)) optimizer%print_every = print_every
     219          408 :       IF (PRESENT(master)) optimizer%master = master
     220          408 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     221              :          !MK This has to be adapted for a new L-BFGS version possibly
     222          204 :          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         1020 :                    optimizer%isave(44))
     225              :          ALLOCATE (optimizer%x(n), optimizer%lower_bound(n), &
     226              :                    optimizer%upper_bound(n), optimizer%gradient(n), &
     227         2244 :                    optimizer%dsave(29), optimizer%work_array(lenwa))
     228        46689 :          optimizer%x = x0
     229          204 :          optimizer%task = 'START'
     230          204 :          optimizer%wanted_relative_f_delta = wanted_relative_f_delta
     231          204 :          optimizer%wanted_projected_gradient = wanted_projected_gradient
     232        46689 :          optimizer%kind_of_bound = 0
     233          204 :          IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound
     234          204 :          IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound
     235          204 :          IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound
     236          204 :          IF (PRESENT(trust_radius)) optimizer%trust_radius = trust_radius
     237              : 
     238              :          CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     239              :                      optimizer%lower_bound, optimizer%upper_bound, &
     240              :                      optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     241              :                      optimizer%wanted_relative_f_delta, &
     242              :                      optimizer%wanted_projected_gradient, optimizer%work_array, &
     243              :                      optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     244              :                      optimizer%csave, optimizer%lsave, optimizer%isave, &
     245          204 :                      optimizer%dsave, optimizer%trust_radius)
     246              :       ELSE
     247              :          NULLIFY ( &
     248          204 :             optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, &
     249          204 :             optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, &
     250          204 :             optimizer%dsave, optimizer%work_array)
     251          612 :          ALLOCATE (optimizer%x(n))
     252        46689 :          optimizer%x(:) = 0.0_dp
     253          612 :          ALLOCATE (optimizer%gradient(n))
     254        46689 :          optimizer%gradient(:) = 0.0_dp
     255              :       END IF
     256       186348 :       CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     257          408 :       optimizer%status = 0
     258              : 
     259          408 :       CALL timestop(handle)
     260              : 
     261          408 :    END SUBROUTINE cp_opt_gopt_create
     262              : 
     263              : ! **************************************************************************************************
     264              : !> \brief releases the optimizer (see doc/ReferenceCounting.html)
     265              : !> \param optimizer the object that should be released
     266              : !> \par History
     267              : !>      02.2002 created [fawzi]
     268              : !>      09.2003 dealloc_ref->release [fawzi]
     269              : !> \author Fawzi Mohamed
     270              : ! **************************************************************************************************
     271          408 :    SUBROUTINE cp_opt_gopt_release(optimizer)
     272              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT)        :: optimizer
     273              : 
     274              :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_release'
     275              : 
     276              :       INTEGER                                            :: handle
     277              : 
     278          408 :       CALL timeset(routineN, handle)
     279              : 
     280          408 :       IF (ASSOCIATED(optimizer%kind_of_bound)) THEN
     281          204 :          DEALLOCATE (optimizer%kind_of_bound)
     282              :       END IF
     283          408 :       IF (ASSOCIATED(optimizer%i_work_array)) THEN
     284          204 :          DEALLOCATE (optimizer%i_work_array)
     285              :       END IF
     286          408 :       IF (ASSOCIATED(optimizer%isave)) THEN
     287          204 :          DEALLOCATE (optimizer%isave)
     288              :       END IF
     289          408 :       IF (ASSOCIATED(optimizer%x)) THEN
     290          408 :          DEALLOCATE (optimizer%x)
     291              :       END IF
     292          408 :       IF (ASSOCIATED(optimizer%lower_bound)) THEN
     293          204 :          DEALLOCATE (optimizer%lower_bound)
     294              :       END IF
     295          408 :       IF (ASSOCIATED(optimizer%upper_bound)) THEN
     296          204 :          DEALLOCATE (optimizer%upper_bound)
     297              :       END IF
     298          408 :       IF (ASSOCIATED(optimizer%gradient)) THEN
     299          408 :          DEALLOCATE (optimizer%gradient)
     300              :       END IF
     301          408 :       IF (ASSOCIATED(optimizer%dsave)) THEN
     302          204 :          DEALLOCATE (optimizer%dsave)
     303              :       END IF
     304          408 :       IF (ASSOCIATED(optimizer%work_array)) THEN
     305          204 :          DEALLOCATE (optimizer%work_array)
     306              :       END IF
     307          408 :       CALL mp_para_env_release(optimizer%para_env)
     308          408 :       CALL gopt_f_release(optimizer%obj_funct)
     309              : 
     310          408 :       CALL timestop(handle)
     311          408 :    END SUBROUTINE cp_opt_gopt_release
     312              : 
     313              : ! **************************************************************************************************
     314              : !> \brief takes different valuse from the optimizer
     315              : !> \param optimizer ...
     316              : !> \param para_env ...
     317              : !> \param obj_funct ...
     318              : !> \param m ...
     319              : !> \param print_every ...
     320              : !> \param wanted_relative_f_delta ...
     321              : !> \param wanted_projected_gradient ...
     322              : !> \param x ...
     323              : !> \param lower_bound ...
     324              : !> \param upper_bound ...
     325              : !> \param kind_of_bound ...
     326              : !> \param master ...
     327              : !> \param actual_projected_gradient ...
     328              : !> \param n_var ...
     329              : !> \param n_iter ...
     330              : !> \param status ...
     331              : !> \param max_f_per_iter ...
     332              : !> \param at_end ...
     333              : !> \param is_master ...
     334              : !> \param last_f ...
     335              : !> \param f ...
     336              : !> \par History
     337              : !>      none
     338              : !> \author Fawzi Mohamed
     339              : !>      @version 2.2002
     340              : ! **************************************************************************************************
     341            0 :    SUBROUTINE cp_opt_gopt_get(optimizer, para_env, &
     342              :                               obj_funct, m, print_every, &
     343              :                               wanted_relative_f_delta, wanted_projected_gradient, &
     344              :                               x, lower_bound, upper_bound, kind_of_bound, master, &
     345              :                               actual_projected_gradient, &
     346              :                               n_var, n_iter, status, max_f_per_iter, at_end, &
     347              :                               is_master, last_f, f)
     348              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN)           :: optimizer
     349              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
     350              :       TYPE(gopt_f_type), OPTIONAL, POINTER               :: obj_funct
     351              :       INTEGER, INTENT(out), OPTIONAL                     :: m, print_every
     352              :       REAL(kind=dp), INTENT(out), OPTIONAL               :: wanted_relative_f_delta, &
     353              :                                                             wanted_projected_gradient
     354              :       REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER     :: x, lower_bound, upper_bound
     355              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: kind_of_bound
     356              :       INTEGER, INTENT(out), OPTIONAL                     :: master
     357              :       REAL(kind=dp), INTENT(out), OPTIONAL               :: actual_projected_gradient
     358              :       INTEGER, INTENT(out), OPTIONAL                     :: n_var, n_iter, status, max_f_per_iter
     359              :       LOGICAL, INTENT(out), OPTIONAL                     :: at_end, is_master
     360              :       REAL(kind=dp), INTENT(out), OPTIONAL               :: last_f, f
     361              : 
     362            0 :       IF (PRESENT(is_master)) is_master = optimizer%master == optimizer%para_env%mepos
     363            0 :       IF (PRESENT(master)) master = optimizer%master
     364            0 :       IF (PRESENT(status)) status = optimizer%status
     365            0 :       IF (PRESENT(para_env)) para_env => optimizer%para_env
     366            0 :       IF (PRESENT(obj_funct)) obj_funct = optimizer%obj_funct
     367            0 :       IF (PRESENT(m)) m = optimizer%m
     368            0 :       IF (PRESENT(max_f_per_iter)) max_f_per_iter = optimizer%max_f_per_iter
     369            0 :       IF (PRESENT(wanted_projected_gradient)) &
     370            0 :          wanted_projected_gradient = optimizer%wanted_projected_gradient
     371            0 :       IF (PRESENT(wanted_relative_f_delta)) &
     372            0 :          wanted_relative_f_delta = optimizer%wanted_relative_f_delta
     373            0 :       IF (PRESENT(print_every)) print_every = optimizer%print_every
     374            0 :       IF (PRESENT(x)) x => optimizer%x
     375            0 :       IF (PRESENT(n_var)) n_var = SIZE(x)
     376            0 :       IF (PRESENT(lower_bound)) lower_bound => optimizer%lower_bound
     377            0 :       IF (PRESENT(upper_bound)) upper_bound => optimizer%upper_bound
     378            0 :       IF (PRESENT(kind_of_bound)) kind_of_bound => optimizer%kind_of_bound
     379            0 :       IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
     380            0 :       IF (PRESENT(last_f)) last_f = optimizer%last_f
     381            0 :       IF (PRESENT(f)) f = optimizer%f
     382            0 :       IF (PRESENT(at_end)) at_end = optimizer%status > 3
     383            0 :       IF (PRESENT(actual_projected_gradient)) &
     384            0 :          actual_projected_gradient = optimizer%projected_gradient
     385            0 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     386            0 :          IF (optimizer%isave(30) > 1 .AND. (optimizer%task(1:5) == "NEW_X" .OR. &
     387              :                                             optimizer%task(1:4) == "STOP" .AND. optimizer%task(7:9) == "CPU")) THEN
     388              :             ! nr iterations >1 .and. dsave contains the wanted data
     389            0 :             IF (PRESENT(last_f)) last_f = optimizer%dsave(2)
     390            0 :             IF (PRESENT(actual_projected_gradient)) &
     391            0 :                actual_projected_gradient = optimizer%dsave(13)
     392              :          ELSE
     393            0 :             CPASSERT(.NOT. PRESENT(last_f))
     394            0 :             CPASSERT(.NOT. PRESENT(actual_projected_gradient))
     395              :          END IF
     396            0 :       ELSE IF (PRESENT(lower_bound) .OR. PRESENT(upper_bound) .OR. PRESENT(kind_of_bound)) THEN
     397            0 :          CPWARN("asked undefined types")
     398              :       END IF
     399              : 
     400            0 :    END SUBROUTINE cp_opt_gopt_get
     401              : 
     402              : ! **************************************************************************************************
     403              : !> \brief does one optimization step
     404              : !> \param optimizer ...
     405              : !> \param n_iter ...
     406              : !> \param f ...
     407              : !> \param last_f ...
     408              : !> \param projected_gradient ...
     409              : !> \param converged ...
     410              : !> \param geo_section ...
     411              : !> \param force_env ...
     412              : !> \param gopt_param ...
     413              : !> \param spgr ...
     414              : !> \par History
     415              : !>      01.2020 modified [pcazade]
     416              : !> \author Fawzi Mohamed
     417              : !>      @version 2.2002
     418              : !> \note
     419              : !>      use directly mainlb in place of setulb ??
     420              : ! **************************************************************************************************
     421         4120 :    SUBROUTINE cp_opt_gopt_step(optimizer, n_iter, f, last_f, &
     422              :                                projected_gradient, converged, geo_section, force_env, &
     423              :                                gopt_param, spgr)
     424              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT)        :: optimizer
     425              :       INTEGER, INTENT(out), OPTIONAL                     :: n_iter
     426              :       REAL(kind=dp), INTENT(out), OPTIONAL               :: f, last_f, projected_gradient
     427              :       LOGICAL, INTENT(out), OPTIONAL                     :: converged
     428              :       TYPE(section_vals_type), POINTER                   :: geo_section
     429              :       TYPE(force_env_type), POINTER                      :: force_env
     430              :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     431              :       TYPE(spgr_type), OPTIONAL, POINTER                 :: spgr
     432              : 
     433              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_opt_gopt_step'
     434              : 
     435              :       CHARACTER(LEN=5)                                   :: wildcard
     436              :       INTEGER                                            :: dataunit, handle, its
     437              :       LOGICAL                                            :: conv, is_master, justEntred, &
     438              :                                                             keep_space_group
     439              :       REAL(KIND=dp)                                      :: t_diff, t_now, t_old
     440         4120 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xold
     441              :       TYPE(cp_logger_type), POINTER                      :: logger
     442              :       TYPE(cp_subsys_type), POINTER                      :: subsys
     443              : 
     444         4120 :       NULLIFY (logger, xold)
     445         8240 :       logger => cp_get_default_logger()
     446         4120 :       CALL timeset(routineN, handle)
     447         4120 :       justEntred = .TRUE.
     448         4120 :       is_master = optimizer%master == optimizer%para_env%mepos
     449         4120 :       IF (PRESENT(converged)) converged = optimizer%status == 4
     450        12360 :       ALLOCATE (xold(SIZE(optimizer%x)))
     451              : 
     452              :       ! collecting subsys
     453         4120 :       CALL force_env_get(force_env, subsys=subsys)
     454              : 
     455         4120 :       keep_space_group = .FALSE.
     456         4120 :       IF (PRESENT(spgr)) THEN
     457         4120 :          IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
     458              :       END IF
     459              : 
     460              :       ! applies rotation matrices to coordinates
     461         4120 :       IF (keep_space_group) THEN
     462            2 :          CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     463              :       END IF
     464              : 
     465      1847524 :       xold = optimizer%x
     466         4120 :       t_old = m_walltime()
     467              : 
     468         4120 :       IF (optimizer%status >= 4) THEN
     469            0 :          CPWARN("status>=4, trying to restart")
     470            0 :          optimizer%status = 0
     471            0 :          IF (is_master) THEN
     472            0 :             optimizer%task = 'START'
     473              :             CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     474              :                         optimizer%lower_bound, optimizer%upper_bound, &
     475              :                         optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     476              :                         optimizer%wanted_relative_f_delta, &
     477              :                         optimizer%wanted_projected_gradient, optimizer%work_array, &
     478              :                         optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     479              :                         optimizer%csave, optimizer%lsave, optimizer%isave, &
     480            0 :                         optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     481              :          END IF
     482              :       END IF
     483              : 
     484              :       DO
     485        13044 :          ifMaster: IF (is_master) THEN
     486         6522 :             IF (optimizer%task(1:7) == 'RESTART') THEN
     487              :                ! restart the optimizer
     488            0 :                optimizer%status = 0
     489            0 :                optimizer%task = 'START'
     490              :                ! applies rotation matrices to coordinates and forces
     491            0 :                IF (keep_space_group) THEN
     492            0 :                   CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     493            0 :                   CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     494              :                END IF
     495              :                CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     496              :                            optimizer%lower_bound, optimizer%upper_bound, &
     497              :                            optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     498              :                            optimizer%wanted_relative_f_delta, &
     499              :                            optimizer%wanted_projected_gradient, optimizer%work_array, &
     500              :                            optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     501              :                            optimizer%csave, optimizer%lsave, optimizer%isave, &
     502            0 :                            optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     503            0 :                IF (keep_space_group) THEN
     504            0 :                   CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     505            0 :                   CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     506              :                END IF
     507              :             END IF
     508         6522 :             IF (optimizer%task(1:2) == 'FG') THEN
     509         2605 :                IF (optimizer%isave(36) > optimizer%max_f_per_iter) THEN
     510            0 :                   optimizer%task = 'STOP: CPU, hit max f eval in iter'
     511            0 :                   optimizer%status = 5 ! anormal exit
     512              :                   CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     513              :                               optimizer%lower_bound, optimizer%upper_bound, &
     514              :                               optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     515              :                               optimizer%wanted_relative_f_delta, &
     516              :                               optimizer%wanted_projected_gradient, optimizer%work_array, &
     517              :                               optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     518              :                               optimizer%csave, optimizer%lsave, optimizer%isave, &
     519            0 :                               optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     520              :                ELSE
     521         2605 :                   optimizer%status = 1
     522              :                END IF
     523         3917 :             ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN
     524         3917 :                IF (justEntred) THEN
     525         1857 :                   optimizer%status = 2
     526              :                   ! applies rotation matrices to coordinates and forces
     527         1857 :                   IF (keep_space_group) THEN
     528            0 :                      CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     529            0 :                      CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     530              :                   END IF
     531              :                   CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     532              :                               optimizer%lower_bound, optimizer%upper_bound, &
     533              :                               optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     534              :                               optimizer%wanted_relative_f_delta, &
     535              :                               optimizer%wanted_projected_gradient, optimizer%work_array, &
     536              :                               optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     537              :                               optimizer%csave, optimizer%lsave, optimizer%isave, &
     538         1857 :                               optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     539         1857 :                   IF (keep_space_group) THEN
     540            0 :                      CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     541            0 :                      CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     542              :                   END IF
     543              :                ELSE
     544              :                   ! applies rotation matrices to coordinates and forces
     545         2060 :                   IF (keep_space_group) THEN
     546            1 :                      CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     547            1 :                      CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     548              :                   END IF
     549         2060 :                   optimizer%status = 3
     550              :                END IF
     551            0 :             ELSE IF (optimizer%task(1:4) == 'CONV') THEN
     552            0 :                optimizer%status = 4
     553            0 :             ELSE IF (optimizer%task(1:4) == 'STOP') THEN
     554            0 :                optimizer%status = 5
     555            0 :                CPWARN("task became stop in an unknown way")
     556            0 :             ELSE IF (optimizer%task(1:5) == 'ERROR') THEN
     557            0 :                optimizer%status = 5
     558              :             ELSE
     559            0 :                CPWARN("unknown task '"//optimizer%task//"'")
     560              :             END IF
     561              :          END IF ifMaster
     562        13044 :          CALL optimizer%para_env%bcast(optimizer%status, optimizer%master)
     563              :          ! Dump info
     564        13044 :          IF (optimizer%status == 3) THEN
     565         4120 :             its = 0
     566         4120 :             IF (is_master) THEN
     567              :                ! Iteration level is taken into account in the optimizer external loop
     568         2060 :                its = optimizer%isave(30)
     569              :             END IF
     570              :          END IF
     571              :          !
     572         5210 :          SELECT CASE (optimizer%status)
     573              :          CASE (1)
     574              :             !op=1 evaluate f and g
     575              :             CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, &
     576              :                             f=optimizer%f, &
     577              :                             gradient=optimizer%gradient, &
     578              :                             final_evaluation=.FALSE., &
     579         5210 :                             master=optimizer%master, para_env=optimizer%para_env)
     580              :             ! do not use keywords?
     581         5210 :             IF (is_master) THEN
     582              :                ! applies rotation matrices to coordinates and forces
     583         2605 :                IF (keep_space_group) THEN
     584            2 :                   CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     585            2 :                   CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     586              :                END IF
     587              :                CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     588              :                            optimizer%lower_bound, optimizer%upper_bound, &
     589              :                            optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     590              :                            optimizer%wanted_relative_f_delta, &
     591              :                            optimizer%wanted_projected_gradient, optimizer%work_array, &
     592              :                            optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     593              :                            optimizer%csave, optimizer%lsave, optimizer%isave, &
     594         2605 :                            optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     595         2605 :                IF (keep_space_group) THEN
     596            2 :                   CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     597            2 :                   CALL spgr_apply_rotations_force(spgr, optimizer%gradient)
     598              :                END IF
     599              :             END IF
     600      4332746 :             CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     601              :          CASE (2)
     602              :             !op=2 begin new iter
     603      3504630 :             CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     604         3714 :             t_old = m_walltime()
     605              :          CASE (3)
     606              :             !op=3 ended iter
     607         4120 :             wildcard = "LBFGS"
     608              :             dataunit = cp_print_key_unit_nr(logger, geo_section, &
     609         4120 :                                             "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
     610         4120 :             IF (is_master) its = optimizer%isave(30)
     611         4120 :             CALL optimizer%para_env%bcast(its, optimizer%master)
     612              : 
     613              :             ! Some IO and Convergence check
     614         4120 :             t_now = m_walltime()
     615         4120 :             t_diff = t_now - t_old
     616         4120 :             t_old = t_now
     617              :             CALL gopt_f_io(optimizer%obj_funct, force_env, force_env%root_section, &
     618              :                            its, optimizer%f, dataunit, optimizer%eold, optimizer%emin, wildcard, gopt_param, &
     619      1847524 :                            SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff)
     620         4120 :             CALL optimizer%para_env%bcast(conv, optimizer%master)
     621              :             CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
     622         4120 :                                               "PRINT%PROGRAM_RUN_INFO")
     623         4120 :             optimizer%eold = optimizer%f
     624         4120 :             optimizer%emin = MIN(optimizer%emin, optimizer%eold)
     625      1847524 :             xold = optimizer%x
     626         4120 :             IF (PRESENT(converged)) converged = conv
     627            0 :             EXIT
     628              :          CASE (4)
     629              :             !op=4 (convergence - normal exit)
     630              :             ! Specific L-BFGS convergence criteria.. overrides the convergence criteria on
     631              :             ! stepsize and gradients
     632              :             dataunit = cp_print_key_unit_nr(logger, geo_section, &
     633            0 :                                             "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
     634            0 :             IF (dataunit > 0) THEN
     635            0 :                WRITE (dataunit, '(T2,A)') ""
     636            0 :                WRITE (dataunit, '(T2,A)') "***********************************************"
     637            0 :                WRITE (dataunit, '(T2,A)') "* Specific L-BFGS convergence criteria         "
     638            0 :                WRITE (dataunit, '(T2,A)') "* WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR  "
     639            0 :                WRITE (dataunit, '(T2,A)') "* satisfied .... run CONVERGED!                "
     640            0 :                WRITE (dataunit, '(T2,A)') "***********************************************"
     641            0 :                WRITE (dataunit, '(T2,A)') ""
     642              :             END IF
     643              :             CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
     644            0 :                                               "PRINT%PROGRAM_RUN_INFO")
     645            0 :             IF (PRESENT(converged)) converged = .TRUE.
     646            0 :             EXIT
     647              :          CASE (5)
     648              :             ! op=5 abnormal exit ()
     649            0 :             CALL optimizer%para_env%bcast(optimizer%task, optimizer%master)
     650              :          CASE (6)
     651              :             ! deallocated
     652            0 :             CPABORT("step on a deallocated opt structure ")
     653              :          CASE default
     654              :             CALL cp_abort(__LOCATION__, &
     655            0 :                           "unknown status "//cp_to_string(optimizer%status))
     656            0 :             optimizer%status = 5
     657        13044 :             EXIT
     658              :          END SELECT
     659         8924 :          IF (optimizer%status == 1 .AND. justEntred) THEN
     660          406 :             optimizer%eold = optimizer%f
     661          406 :             optimizer%emin = optimizer%eold
     662              :          END IF
     663              :          justEntred = .FALSE.
     664              :       END DO
     665              : 
     666      3690928 :       CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     667              :       CALL cp_opt_gopt_bcast_res(optimizer, &
     668              :                                  n_iter=optimizer%n_iter, &
     669              :                                  f=optimizer%f, last_f=optimizer%last_f, &
     670         4120 :                                  projected_gradient=optimizer%projected_gradient)
     671              : 
     672         4120 :       DEALLOCATE (xold)
     673         4120 :       IF (PRESENT(f)) f = optimizer%f
     674         4120 :       IF (PRESENT(last_f)) last_f = optimizer%last_f
     675         4120 :       IF (PRESENT(projected_gradient)) &
     676            0 :          projected_gradient = optimizer%projected_gradient
     677         4120 :       IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
     678         4120 :       CALL timestop(handle)
     679              : 
     680         4120 :    END SUBROUTINE cp_opt_gopt_step
     681              : 
     682              : ! **************************************************************************************************
     683              : !> \brief returns the results (and broadcasts them)
     684              : !> \param optimizer the optimizer object the info is taken from
     685              : !> \param n_iter the number of iterations
     686              : !> \param f the actual value of the objective function (f)
     687              : !> \param last_f the last value of f
     688              : !> \param projected_gradient the infinity norm of the projected gradient
     689              : !> \par History
     690              : !>      none
     691              : !> \author Fawzi Mohamed
     692              : !>      @version 2.2002
     693              : !> \note
     694              : !>      private routine
     695              : ! **************************************************************************************************
     696         4120 :    SUBROUTINE cp_opt_gopt_bcast_res(optimizer, n_iter, f, last_f, &
     697              :                                     projected_gradient)
     698              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(IN)           :: optimizer
     699              :       INTEGER, INTENT(out), OPTIONAL                     :: n_iter
     700              :       REAL(kind=dp), INTENT(inout), OPTIONAL             :: f, last_f, projected_gradient
     701              : 
     702              :       REAL(kind=dp), DIMENSION(4)                        :: results
     703              : 
     704         4120 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     705              :          results = (/REAL(optimizer%isave(30), kind=dp), &
     706        10300 :                      optimizer%f, optimizer%dsave(2), optimizer%dsave(13)/)
     707              :       END IF
     708         4120 :       CALL optimizer%para_env%bcast(results, optimizer%master)
     709         4120 :       IF (PRESENT(n_iter)) n_iter = NINT(results(1))
     710         4120 :       IF (PRESENT(f)) f = results(2)
     711         4120 :       IF (PRESENT(last_f)) last_f = results(3)
     712         4120 :       IF (PRESENT(projected_gradient)) &
     713         4120 :          projected_gradient = results(4)
     714              : 
     715         4120 :    END SUBROUTINE cp_opt_gopt_bcast_res
     716              : 
     717              : ! **************************************************************************************************
     718              : !> \brief goes to the next optimal point (after an optimizer iteration)
     719              : !>      returns true if converged
     720              : !> \param optimizer the optimizer that goes to the next point
     721              : !> \param n_iter ...
     722              : !> \param f ...
     723              : !> \param last_f ...
     724              : !> \param projected_gradient ...
     725              : !> \param converged ...
     726              : !> \param geo_section ...
     727              : !> \param force_env ...
     728              : !> \param gopt_param ...
     729              : !> \param spgr ...
     730              : !> \return ...
     731              : !> \par History
     732              : !>      01.2020 modified [pcazade]
     733              : !> \author Fawzi Mohamed
     734              : !>      @version 2.2002
     735              : !> \note
     736              : !>      if you deactivate convergence control it returns never false
     737              : ! **************************************************************************************************
     738         4120 :    FUNCTION cp_opt_gopt_next(optimizer, n_iter, f, last_f, &
     739              :                              projected_gradient, converged, geo_section, force_env, &
     740              :                              gopt_param, spgr) RESULT(res)
     741              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT)        :: optimizer
     742              :       INTEGER, INTENT(out), OPTIONAL                     :: n_iter
     743              :       REAL(kind=dp), INTENT(out), OPTIONAL               :: f, last_f, projected_gradient
     744              :       LOGICAL, INTENT(out)                               :: converged
     745              :       TYPE(section_vals_type), POINTER                   :: geo_section
     746              :       TYPE(force_env_type), POINTER                      :: force_env
     747              :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     748              :       TYPE(spgr_type), OPTIONAL, POINTER                 :: spgr
     749              :       LOGICAL                                            :: res
     750              : 
     751              :       ! passes spgr structure if present
     752              :       CALL cp_opt_gopt_step(optimizer, n_iter=n_iter, f=f, &
     753              :                             last_f=last_f, projected_gradient=projected_gradient, &
     754              :                             converged=converged, geo_section=geo_section, &
     755         4120 :                             force_env=force_env, gopt_param=gopt_param, spgr=spgr)
     756         4120 :       res = (optimizer%status < 40) .AND. .NOT. converged
     757              : 
     758         4120 :    END FUNCTION cp_opt_gopt_next
     759              : 
     760              : ! **************************************************************************************************
     761              : !> \brief stops the optimization
     762              : !> \param optimizer ...
     763              : !> \par History
     764              : !>      none
     765              : !> \author Fawzi Mohamed
     766              : !>      @version 2.2002
     767              : ! **************************************************************************************************
     768            0 :    SUBROUTINE cp_opt_gopt_stop(optimizer)
     769              :       TYPE(cp_lbfgs_opt_gopt_type), INTENT(INOUT)       :: optimizer
     770              : 
     771            0 :       optimizer%task = 'STOPPED on user request'
     772            0 :       optimizer%status = 4 ! normal exit
     773            0 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     774              :          CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, &
     775              :                      optimizer%lower_bound, optimizer%upper_bound, &
     776              :                      optimizer%kind_of_bound, optimizer%f, optimizer%gradient, &
     777              :                      optimizer%wanted_relative_f_delta, &
     778              :                      optimizer%wanted_projected_gradient, optimizer%work_array, &
     779              :                      optimizer%i_work_array, optimizer%task, optimizer%print_every, &
     780              :                      optimizer%csave, optimizer%lsave, optimizer%isave, &
     781            0 :                      optimizer%dsave, optimizer%trust_radius)
     782              :       END IF
     783              : 
     784            0 :    END SUBROUTINE cp_opt_gopt_stop
     785              : 
     786            0 : END MODULE cp_lbfgs_optimizer_gopt
        

Generated by: LCOV version 2.0-1