LCOV - code coverage report
Current view: top level - src/motion - cp_lbfgs_optimizer_gopt.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 169 259 65.3 %
Date: 2024-12-12 07:20:04 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        2352 :    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         392 :       CALL timeset(routineN, handle)
     191             : 
     192             :       NULLIFY (optimizer%kind_of_bound, &
     193         392 :                optimizer%i_work_array, &
     194         392 :                optimizer%isave, &
     195         392 :                optimizer%x, &
     196         392 :                optimizer%lower_bound, &
     197         392 :                optimizer%upper_bound, &
     198         392 :                optimizer%gradient, &
     199         392 :                optimizer%dsave, &
     200         392 :                optimizer%work_array, &
     201             :                optimizer%para_env, &
     202         392 :                optimizer%obj_funct)
     203         392 :       n = SIZE(x0)
     204         392 :       optimizer%m = 4
     205         392 :       IF (PRESENT(m)) optimizer%m = m
     206         392 :       optimizer%master = para_env%source
     207         392 :       optimizer%para_env => para_env
     208         392 :       CALL para_env%retain()
     209         392 :       optimizer%obj_funct => obj_funct
     210         392 :       CALL gopt_f_retain(obj_funct)
     211         392 :       optimizer%max_f_per_iter = 20
     212         392 :       IF (PRESENT(max_f_per_iter)) optimizer%max_f_per_iter = max_f_per_iter
     213         392 :       optimizer%print_every = -1
     214         392 :       optimizer%n_iter = 0
     215         392 :       optimizer%f = -1.0_dp
     216         392 :       optimizer%last_f = -1.0_dp
     217         392 :       optimizer%projected_gradient = -1.0_dp
     218         392 :       IF (PRESENT(print_every)) optimizer%print_every = print_every
     219         392 :       IF (PRESENT(master)) optimizer%master = master
     220         392 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     221             :          !MK This has to be adapted for a new L-BFGS version possibly
     222         196 :          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         980 :                    optimizer%isave(44))
     225             :          ALLOCATE (optimizer%x(n), optimizer%lower_bound(n), &
     226             :                    optimizer%upper_bound(n), optimizer%gradient(n), &
     227        2156 :                    optimizer%dsave(29), optimizer%work_array(lenwa))
     228       45709 :          optimizer%x = x0
     229         196 :          optimizer%task = 'START'
     230         196 :          optimizer%wanted_relative_f_delta = wanted_relative_f_delta
     231         196 :          optimizer%wanted_projected_gradient = wanted_projected_gradient
     232       45709 :          optimizer%kind_of_bound = 0
     233         196 :          IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound
     234         196 :          IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound
     235         196 :          IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound
     236         196 :          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         196 :                      optimizer%dsave, optimizer%trust_radius)
     246             :       ELSE
     247             :          NULLIFY ( &
     248         196 :             optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, &
     249         196 :             optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, &
     250         196 :             optimizer%dsave, optimizer%work_array)
     251         588 :          ALLOCATE (optimizer%x(n))
     252       45709 :          optimizer%x(:) = 0.0_dp
     253         588 :          ALLOCATE (optimizer%gradient(n))
     254       45709 :          optimizer%gradient(:) = 0.0_dp
     255             :       END IF
     256      182444 :       CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     257         392 :       optimizer%status = 0
     258             : 
     259         392 :       CALL timestop(handle)
     260             : 
     261         392 :    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         392 :    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         392 :       CALL timeset(routineN, handle)
     279             : 
     280         392 :       IF (ASSOCIATED(optimizer%kind_of_bound)) THEN
     281         196 :          DEALLOCATE (optimizer%kind_of_bound)
     282             :       END IF
     283         392 :       IF (ASSOCIATED(optimizer%i_work_array)) THEN
     284         196 :          DEALLOCATE (optimizer%i_work_array)
     285             :       END IF
     286         392 :       IF (ASSOCIATED(optimizer%isave)) THEN
     287         196 :          DEALLOCATE (optimizer%isave)
     288             :       END IF
     289         392 :       IF (ASSOCIATED(optimizer%x)) THEN
     290         392 :          DEALLOCATE (optimizer%x)
     291             :       END IF
     292         392 :       IF (ASSOCIATED(optimizer%lower_bound)) THEN
     293         196 :          DEALLOCATE (optimizer%lower_bound)
     294             :       END IF
     295         392 :       IF (ASSOCIATED(optimizer%upper_bound)) THEN
     296         196 :          DEALLOCATE (optimizer%upper_bound)
     297             :       END IF
     298         392 :       IF (ASSOCIATED(optimizer%gradient)) THEN
     299         392 :          DEALLOCATE (optimizer%gradient)
     300             :       END IF
     301         392 :       IF (ASSOCIATED(optimizer%dsave)) THEN
     302         196 :          DEALLOCATE (optimizer%dsave)
     303             :       END IF
     304         392 :       IF (ASSOCIATED(optimizer%work_array)) THEN
     305         196 :          DEALLOCATE (optimizer%work_array)
     306             :       END IF
     307         392 :       CALL mp_para_env_release(optimizer%para_env)
     308         392 :       CALL gopt_f_release(optimizer%obj_funct)
     309             : 
     310         392 :       CALL timestop(handle)
     311         392 :    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        4110 :    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        4110 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xold
     441             :       TYPE(cp_logger_type), POINTER                      :: logger
     442             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     443             : 
     444        4110 :       NULLIFY (logger, xold)
     445        8220 :       logger => cp_get_default_logger()
     446        4110 :       CALL timeset(routineN, handle)
     447        4110 :       justEntred = .TRUE.
     448        4110 :       is_master = optimizer%master == optimizer%para_env%mepos
     449        4110 :       IF (PRESENT(converged)) converged = optimizer%status == 4
     450       12330 :       ALLOCATE (xold(SIZE(optimizer%x)))
     451             : 
     452             :       ! collecting subsys
     453        4110 :       CALL force_env_get(force_env, subsys=subsys)
     454             : 
     455        4110 :       keep_space_group = .FALSE.
     456        4110 :       IF (PRESENT(spgr)) THEN
     457        4110 :          IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
     458             :       END IF
     459             : 
     460             :       ! applies rotation matrices to coordinates
     461        4110 :       IF (keep_space_group) THEN
     462           2 :          CALL spgr_apply_rotations_coord(spgr, optimizer%x)
     463             :       END IF
     464             : 
     465     1849302 :       xold = optimizer%x
     466        4110 :       t_old = m_walltime()
     467             : 
     468        4110 :       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       12982 :          ifMaster: IF (is_master) THEN
     486        6491 :             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        6491 :             IF (optimizer%task(1:2) == 'FG') THEN
     509        2576 :                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        2576 :                   optimizer%status = 1
     522             :                END IF
     523        3915 :             ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN
     524        3915 :                IF (justEntred) THEN
     525        1860 :                   optimizer%status = 2
     526             :                   ! applies rotation matrices to coordinates and forces
     527        1860 :                   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        1860 :                               optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     539        1860 :                   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        2055 :                   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        2055 :                   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       12982 :          CALL optimizer%para_env%bcast(optimizer%status, optimizer%master)
     563             :          ! Dump info
     564       12982 :          IF (optimizer%status == 3) THEN
     565        4110 :             its = 0
     566        4110 :             IF (is_master) THEN
     567             :                ! Iteration level is taken into account in the optimizer external loop
     568        2055 :                its = optimizer%isave(30)
     569             :             END IF
     570             :          END IF
     571             :          !
     572        5152 :          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        5152 :                             master=optimizer%master, para_env=optimizer%para_env)
     580             :             ! do not use keywords?
     581        5152 :             IF (is_master) THEN
     582             :                ! applies rotation matrices to coordinates and forces
     583        2576 :                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        2576 :                            optimizer%dsave, optimizer%trust_radius, spgr=spgr)
     595        2576 :                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     4256008 :             CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     601             :          CASE (2)
     602             :             !op=2 begin new iter
     603     3512100 :             CALL optimizer%para_env%bcast(optimizer%x, optimizer%master)
     604        3720 :             t_old = m_walltime()
     605             :          CASE (3)
     606             :             !op=3 ended iter
     607        4110 :             wildcard = "LBFGS"
     608             :             dataunit = cp_print_key_unit_nr(logger, geo_section, &
     609        4110 :                                             "PRINT%PROGRAM_RUN_INFO", extension=".geoLog")
     610        4110 :             IF (is_master) its = optimizer%isave(30)
     611        4110 :             CALL optimizer%para_env%bcast(its, optimizer%master)
     612             : 
     613             :             ! Some IO and Convergence check
     614        4110 :             t_now = m_walltime()
     615        4110 :             t_diff = t_now - t_old
     616        4110 :             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     1849302 :                            SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff)
     620        4110 :             CALL optimizer%para_env%bcast(conv, optimizer%master)
     621             :             CALL cp_print_key_finished_output(dataunit, logger, geo_section, &
     622        4110 :                                               "PRINT%PROGRAM_RUN_INFO")
     623        4110 :             optimizer%eold = optimizer%f
     624        4110 :             optimizer%emin = MIN(optimizer%emin, optimizer%eold)
     625     1849302 :             xold = optimizer%x
     626        4110 :             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       12982 :             EXIT
     658             :          END SELECT
     659        8872 :          IF (optimizer%status == 1 .AND. justEntred) THEN
     660         390 :             optimizer%eold = optimizer%f
     661         390 :             optimizer%emin = optimizer%eold
     662             :          END IF
     663             :          justEntred = .FALSE.
     664             :       END DO
     665             : 
     666     3694494 :       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        4110 :                                  projected_gradient=optimizer%projected_gradient)
     671             : 
     672        4110 :       DEALLOCATE (xold)
     673        4110 :       IF (PRESENT(f)) f = optimizer%f
     674        4110 :       IF (PRESENT(last_f)) last_f = optimizer%last_f
     675        4110 :       IF (PRESENT(projected_gradient)) &
     676           0 :          projected_gradient = optimizer%projected_gradient
     677        4110 :       IF (PRESENT(n_iter)) n_iter = optimizer%n_iter
     678        4110 :       CALL timestop(handle)
     679             : 
     680        4110 :    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        4110 :    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        4110 :       IF (optimizer%master == optimizer%para_env%mepos) THEN
     705             :          results = (/REAL(optimizer%isave(30), kind=dp), &
     706       10275 :                      optimizer%f, optimizer%dsave(2), optimizer%dsave(13)/)
     707             :       END IF
     708        4110 :       CALL optimizer%para_env%bcast(results, optimizer%master)
     709        4110 :       IF (PRESENT(n_iter)) n_iter = NINT(results(1))
     710        4110 :       IF (PRESENT(f)) f = results(2)
     711        4110 :       IF (PRESENT(last_f)) last_f = results(3)
     712        4110 :       IF (PRESENT(projected_gradient)) &
     713        4110 :          projected_gradient = results(4)
     714             : 
     715        4110 :    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        4110 :    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        4110 :                             force_env=force_env, gopt_param=gopt_param, spgr=spgr)
     756        4110 :       res = (optimizer%status < 40) .AND. .NOT. converged
     757             : 
     758        4110 :    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 1.15