LCOV - code coverage report
Current view: top level - src/motion - cp_lbfgs_optimizer_gopt.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 169 259 65.3 %
Date: 2024-04-20 06:29:22 Functions: 5 8 62.5 %

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

Generated by: LCOV version 1.15