LCOV - code coverage report
Current view: top level - src/motion - cg_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 458 495 92.5 %
Date: 2023-03-30 11:55:16 Functions: 12 12 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utilities for Geometry optimization using  Conjugate Gradients
      10             : !> \author Teodoro Laino [teo]
      11             : !>      10.2005
      12             : ! **************************************************************************************************
      13             : MODULE cg_utils
      14             :    USE cp_external_control, ONLY: external_control
      15             :    USE dimer_types, ONLY: dimer_env_type
      16             :    USE dimer_utils, ONLY: dimer_thrs, &
      17             :                           rotate_dimer
      18             :    USE global_types, ONLY: global_environment_type
      19             :    USE gopt_f_types, ONLY: gopt_f_type
      20             :    USE gopt_param_types, ONLY: gopt_param_type
      21             :    USE input_constants, ONLY: default_cell_method_id, &
      22             :                               default_minimization_method_id, &
      23             :                               default_shellcore_method_id, &
      24             :                               default_ts_method_id, &
      25             :                               ls_2pnt, &
      26             :                               ls_fit, &
      27             :                               ls_gold
      28             :    USE kinds, ONLY: dp
      29             :    USE mathconstants, ONLY: pi
      30             :    USE memory_utilities, ONLY: reallocate
      31             :    USE message_passing, ONLY: mp_para_env_type
      32             : #include "../base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             :    PRIVATE
      36             : 
      37             :    #:include "gopt_f77_methods.fypp"
      38             : 
      39             :    PUBLIC :: cg_linmin, get_conjugate_direction
      40             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      41             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_utils'
      42             : 
      43             : CONTAINS
      44             : 
      45             : ! **************************************************************************************************
      46             : !> \brief Main driver for line minimization routines for CG
      47             : !> \param gopt_env ...
      48             : !> \param xvec ...
      49             : !> \param xi ...
      50             : !> \param g ...
      51             : !> \param opt_energy ...
      52             : !> \param output_unit ...
      53             : !> \param gopt_param ...
      54             : !> \param globenv ...
      55             : !> \par History
      56             : !>      10.2005 created [tlaino]
      57             : !> \author Teodoro Laino
      58             : ! **************************************************************************************************
      59        2532 :    RECURSIVE SUBROUTINE cg_linmin(gopt_env, xvec, xi, g, opt_energy, output_unit, gopt_param, &
      60             :                                   globenv)
      61             : 
      62             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
      63             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xvec, xi, g
      64             :       REAL(KIND=dp), INTENT(INOUT)                       :: opt_energy
      65             :       INTEGER                                            :: output_unit
      66             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
      67             :       TYPE(global_environment_type), POINTER             :: globenv
      68             : 
      69             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cg_linmin'
      70             : 
      71             :       INTEGER                                            :: handle
      72             : 
      73        2532 :       CALL timeset(routineN, handle)
      74        2532 :       gopt_env%do_line_search = .TRUE.
      75        3446 :       SELECT CASE (gopt_env%type_id)
      76             :       CASE (default_minimization_method_id)
      77        2556 :          SELECT CASE (gopt_param%cg_ls%type_id)
      78             :          CASE (ls_2pnt)
      79             :             CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=gopt_param%cg_ls%grad_only, &
      80         788 :                              output_unit=output_unit)
      81             :          CASE (ls_fit)
      82             :             CALL linmin_fit(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brack_limit, &
      83          54 :                             gopt_param%cg_ls%initial_step, output_unit, gopt_param, globenv)
      84             :          CASE (ls_gold)
      85             :             CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
      86             :                              gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
      87          72 :                              gopt_param%cg_ls%initial_step, output_unit, globenv)
      88             :          CASE DEFAULT
      89         914 :             CPABORT("Line Search type not yet implemented in CG.")
      90             :          END SELECT
      91             :       CASE (default_ts_method_id)
      92        2302 :          SELECT CASE (gopt_param%cg_ls%type_id)
      93             :          CASE (ls_2pnt)
      94         854 :             IF (gopt_env%dimer_rotation) THEN
      95         726 :                CALL rotmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy)
      96             :             ELSE
      97             :                CALL tslmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy, gopt_param, &
      98         128 :                                 output_unit)
      99             :             END IF
     100             :          CASE DEFAULT
     101         854 :             CPABORT("Line Search type not yet implemented in CG for TS search.")
     102             :          END SELECT
     103             :       CASE (default_cell_method_id)
     104        1248 :          SELECT CASE (gopt_param%cg_ls%type_id)
     105             :          CASE (ls_2pnt)
     106             :             CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.TRUE., &
     107         484 :                              output_unit=output_unit)
     108             :          CASE (ls_gold)
     109             :             CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
     110             :                              gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
     111         110 :                              gopt_param%cg_ls%initial_step, output_unit, globenv)
     112             :          CASE DEFAULT
     113         594 :             CPABORT("Line Search type not yet implemented in CG for cell optimization.")
     114             :          END SELECT
     115             :       CASE (default_shellcore_method_id)
     116        2532 :          SELECT CASE (gopt_param%cg_ls%type_id)
     117             :          CASE (ls_2pnt)
     118             :             CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.TRUE., &
     119         170 :                              output_unit=output_unit)
     120             :          CASE DEFAULT
     121         170 :             CPABORT("Line Search type not yet implemented in CG for shellcore optimization.")
     122             :          END SELECT
     123             : 
     124             :       END SELECT
     125        2532 :       gopt_env%do_line_search = .FALSE.
     126        2532 :       CALL timestop(handle)
     127             : 
     128        2532 :    END SUBROUTINE cg_linmin
     129             : 
     130             : ! **************************************************************************************************
     131             : !> \brief Line search subroutine based on 2 points (using gradients and energies
     132             : !>        or only gradients)
     133             : !> \param gopt_env ...
     134             : !> \param x0 ...
     135             : !> \param ls_vec ...
     136             : !> \param g ...
     137             : !> \param opt_energy ...
     138             : !> \param gopt_param ...
     139             : !> \param use_only_grad ...
     140             : !> \param output_unit ...
     141             : !> \author Teodoro Laino - created [tlaino] - 03.2008
     142             : ! **************************************************************************************************
     143        1442 :    RECURSIVE SUBROUTINE linmin_2pnt(gopt_env, x0, ls_vec, g, opt_energy, gopt_param, use_only_grad, &
     144             :                                     output_unit)
     145             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     146             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0, ls_vec, g
     147             :       REAL(KIND=dp), INTENT(INOUT)                       :: opt_energy
     148             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     149             :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_only_grad
     150             :       INTEGER, INTENT(IN)                                :: output_unit
     151             : 
     152             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'linmin_2pnt'
     153             : 
     154             :       INTEGER                                            :: handle
     155             :       LOGICAL                                            :: my_use_only_grad, &
     156             :                                                             save_consistent_energy_force
     157             :       REAL(KIND=dp)                                      :: a, b, c, dx, dx_min, dx_min_save, &
     158             :                                                             dx_thrs, norm_grad1, norm_grad2, &
     159             :                                                             norm_ls_vec, opt_energy2, x_grad_zero
     160        1442 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient2, ls_norm
     161             : 
     162        1442 :       CALL timeset(routineN, handle)
     163      363926 :       norm_ls_vec = SQRT(DOT_PRODUCT(ls_vec, ls_vec))
     164        1442 :       my_use_only_grad = .FALSE.
     165        1442 :       IF (PRESENT(use_only_grad)) my_use_only_grad = use_only_grad
     166        1442 :       IF (norm_ls_vec /= 0.0_dp) THEN
     167        4326 :          ALLOCATE (ls_norm(SIZE(ls_vec)))
     168        4326 :          ALLOCATE (gradient2(SIZE(ls_vec)))
     169      726410 :          ls_norm = ls_vec/norm_ls_vec
     170        1442 :          dx = norm_ls_vec
     171        1442 :          dx_thrs = gopt_param%cg_ls%max_step
     172             : 
     173      726410 :          x0 = x0 + dx*ls_norm
     174             :          ![NB] don't need consistent energies and forces if using only gradient
     175        1442 :          save_consistent_energy_force = gopt_env%require_consistent_energy_force
     176        1442 :          gopt_env%require_consistent_energy_force = .NOT. my_use_only_grad
     177             :          CALL cp_eval_at(gopt_env, x0, opt_energy2, gradient2, master=gopt_env%force_env%para_env%mepos, &
     178        1442 :                          final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     179        1442 :          gopt_env%require_consistent_energy_force = save_consistent_energy_force
     180             : 
     181      363926 :          norm_grad1 = -DOT_PRODUCT(g, ls_norm)
     182      363926 :          norm_grad2 = DOT_PRODUCT(gradient2, ls_norm)
     183        1442 :          IF (my_use_only_grad) THEN
     184             :             ! a*x+b=y
     185             :             ! per x=0; b=norm_grad1
     186         654 :             b = norm_grad1
     187             :             ! per x=dx; a*dx+b=norm_grad2
     188         654 :             a = (norm_grad2 - b)/dx
     189         654 :             x_grad_zero = -b/a
     190         654 :             dx_min = x_grad_zero
     191             :          ELSE
     192             :             ! ax**2+b*x+c=y
     193             :             ! per x=0 ; c=opt_energy
     194         788 :             c = opt_energy
     195             :             ! per x=dx;          a*dx**2 + b*dx + c = opt_energy2
     196             :             ! per x=dx;        2*a*dx    + b        = norm_grad2
     197             :             !
     198             :             !                  - a*dx**2        + c = (opt_energy2-norm_grad2*dx)
     199             :             !                    a*dx**2            = c - (opt_energy2-norm_grad2*dx)
     200         788 :             a = (c - (opt_energy2 - norm_grad2*dx))/dx**2
     201         788 :             b = norm_grad2 - 2.0_dp*a*dx
     202         788 :             dx_min = 0.0_dp
     203         788 :             IF (a /= 0.0_dp) dx_min = -b/(2.0_dp*a)
     204         788 :             opt_energy = opt_energy2
     205             :          END IF
     206        1442 :          dx_min_save = dx_min
     207             :          ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
     208             :          ! step length
     209        1442 :          IF (ABS(dx_min) > dx_thrs) dx_min = SIGN(1.0_dp, dx_min)*dx_thrs
     210      726410 :          x0 = x0 + (dx_min - dx)*ls_norm
     211             : 
     212             :          ! Print out LS info
     213        1442 :          IF (output_unit > 0) THEN
     214         721 :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     215             :             WRITE (UNIT=output_unit, FMT="(T2,A,T31,A,T78,A)") &
     216         721 :                "***", "2PNT LINE SEARCH INFO", "***"
     217         721 :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A)") "***", "***"
     218             :             WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
     219         721 :                "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
     220             :             WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
     221         721 :                "***", "DX (FITTED   )=", dx_min_save, "DX (ACCEPTED )=", dx_min, "***"
     222         721 :             WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     223             :          END IF
     224        1442 :          DEALLOCATE (ls_norm)
     225        1442 :          DEALLOCATE (gradient2)
     226             :       ELSE
     227             :          ! Do Nothing, since.. if the effective force is 0 means that we are already
     228             :          ! in the saddle point..
     229             :       END IF
     230        1442 :       CALL timestop(handle)
     231        1442 :    END SUBROUTINE linmin_2pnt
     232             : 
     233             : ! **************************************************************************************************
     234             : !> \brief Translational minimization for the Dimer Method - 2pnt LS
     235             : !> \param gopt_env ...
     236             : !> \param dimer_env ...
     237             : !> \param x0 ...
     238             : !> \param tls_vec ...
     239             : !> \param opt_energy ...
     240             : !> \param gopt_param ...
     241             : !> \param output_unit ...
     242             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     243             : ! **************************************************************************************************
     244         128 :    SUBROUTINE tslmin_2pnt(gopt_env, dimer_env, x0, tls_vec, opt_energy, gopt_param, output_unit)
     245             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     246             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     247             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0, tls_vec
     248             :       REAL(KIND=dp), INTENT(INOUT)                       :: opt_energy
     249             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     250             :       INTEGER, INTENT(IN)                                :: output_unit
     251             : 
     252             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'tslmin_2pnt'
     253             : 
     254             :       INTEGER                                            :: handle
     255             :       REAL(KIND=dp)                                      :: dx, dx_min, dx_min_acc, dx_min_save, &
     256             :                                                             dx_thrs, norm_tls_vec, opt_energy2
     257         128 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tls_norm
     258             : 
     259         128 :       CALL timeset(routineN, handle)
     260        1886 :       norm_tls_vec = SQRT(DOT_PRODUCT(tls_vec, tls_vec))
     261         128 :       IF (norm_tls_vec /= 0.0_dp) THEN
     262         384 :          ALLOCATE (tls_norm(SIZE(tls_vec)))
     263             : 
     264        3644 :          tls_norm = tls_vec/norm_tls_vec
     265         128 :          dimer_env%tsl%tls_vec => tls_norm
     266             : 
     267         128 :          dx = norm_tls_vec
     268         128 :          dx_thrs = gopt_param%cg_ls%max_step
     269             :          ! If curvature is positive let's make the largest step allowed
     270         128 :          IF (dimer_env%rot%curvature > 0) dx = dx_thrs
     271        3644 :          x0 = x0 + dx*tls_norm
     272             :          CALL cp_eval_at(gopt_env, x0, opt_energy2, master=gopt_env%force_env%para_env%mepos, &
     273         128 :                          final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     274         128 :          IF (dimer_env%rot%curvature > 0) THEN
     275          30 :             dx_min = 0.0_dp
     276          30 :             dx_min_save = dx
     277          30 :             dx_min_acc = dx
     278             :          ELSE
     279             :             ! First let's try to interpolate the minimum
     280          98 :             dx_min = -opt_energy/(opt_energy2 - opt_energy)*dx
     281             :             ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
     282             :             ! step length
     283          98 :             dx_min_save = dx_min
     284          98 :             IF (ABS(dx_min) > dx_thrs) dx_min = SIGN(1.0_dp, dx_min)*dx_thrs
     285          98 :             dx_min_acc = dx_min
     286          98 :             dx_min = dx_min - dx
     287             :          END IF
     288        3644 :          x0 = x0 + dx_min*tls_norm
     289             : 
     290             :          ! Print out LS info
     291         128 :          IF (output_unit > 0) THEN
     292          64 :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     293             :             WRITE (UNIT=output_unit, FMT="(T2,A,T24,A,T78,A)") &
     294          64 :                "***", "2PNT TRANSLATIONAL LINE SEARCH INFO", "***"
     295          64 :             WRITE (UNIT=output_unit, FMT="(T2,A,T78,A)") "***", "***"
     296             :             WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T78,A)") &
     297          64 :                "***", "LOCAL CURVATURE =", dimer_env%rot%curvature, "***"
     298             :             WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
     299          64 :                "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
     300             :             WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
     301          64 :                "***", "DX (FITTED   )=", dx_min_save, "DX (ACCEPTED )=", dx_min_acc, "***"
     302          64 :             WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     303             :          END IF
     304             : 
     305             :          ! Here we compute the value of the energy in point zero..
     306             :          CALL cp_eval_at(gopt_env, x0, opt_energy, master=gopt_env%force_env%para_env%mepos, &
     307         128 :                          final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     308             : 
     309         128 :          DEALLOCATE (tls_norm)
     310             :       ELSE
     311             :          ! Do Nothing, since.. if the effective force is 0 means that we are already
     312             :          ! in the saddle point..
     313             :       END IF
     314         128 :       CALL timestop(handle)
     315             : 
     316         128 :    END SUBROUTINE tslmin_2pnt
     317             : 
     318             : ! **************************************************************************************************
     319             : !> \brief Rotational minimization for the Dimer Method - 2 pnt LS
     320             : !> \param gopt_env ...
     321             : !> \param dimer_env ...
     322             : !> \param x0 ...
     323             : !> \param theta ...
     324             : !> \param opt_energy ...
     325             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     326             : ! **************************************************************************************************
     327         726 :    SUBROUTINE rotmin_2pnt(gopt_env, dimer_env, x0, theta, opt_energy)
     328             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     329             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     330             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0, theta
     331             :       REAL(KIND=dp), INTENT(INOUT)                       :: opt_energy
     332             : 
     333             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rotmin_2pnt'
     334             : 
     335             :       INTEGER                                            :: handle
     336             :       REAL(KIND=dp)                                      :: a0, a1, angle, b1, curvature0, &
     337             :                                                             curvature1, curvature2, dCdp, f
     338         726 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: work
     339             : 
     340         726 :       CALL timeset(routineN, handle)
     341         726 :       curvature0 = dimer_env%rot%curvature
     342         726 :       dCdp = dimer_env%rot%dCdp
     343         726 :       b1 = 0.5_dp*dCdp
     344         726 :       angle = -0.5_dp*ATAN(dCdp/(2.0_dp*ABS(curvature0)))
     345         726 :       dimer_env%rot%angle1 = angle
     346       12840 :       dimer_env%cg_rot%nvec_old = dimer_env%nvec
     347         726 :       IF (angle > dimer_env%rot%angle_tol) THEN
     348             :          ! Rotating the dimer of dtheta degrees
     349         702 :          CALL rotate_dimer(dimer_env%nvec, theta, angle)
     350             :          ! Re-compute energy, gradients and rotation vector for new R1
     351             :          CALL cp_eval_at(gopt_env, x0, f, master=gopt_env%force_env%para_env%mepos, &
     352         702 :                          final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     353             : 
     354         702 :          curvature1 = dimer_env%rot%curvature
     355         702 :          a1 = (curvature0 - curvature1 + b1*SIN(2.0_dp*angle))/(1.0_dp - COS(2.0_dp*angle))
     356         702 :          a0 = 2.0_dp*(curvature0 - a1)
     357         702 :          angle = 0.5_dp*ATAN(b1/a1)
     358         702 :          curvature2 = a0/2.0_dp + a1*COS(2.0_dp*angle) + b1*SIN(2.0_dp*angle)
     359         702 :          IF (curvature2 > curvature0) THEN
     360           4 :             angle = angle + pi/2.0_dp
     361           4 :             curvature2 = a0/2.0_dp + a1*COS(2.0_dp*angle) + b1*SIN(2.0_dp*angle)
     362             :          END IF
     363         702 :          dimer_env%rot%angle2 = angle
     364         702 :          dimer_env%rot%curvature = curvature2
     365             :          ! Rotating the dimer the optimized (in plane) vector position
     366       12708 :          dimer_env%nvec = dimer_env%cg_rot%nvec_old
     367         702 :          CALL rotate_dimer(dimer_env%nvec, theta, angle)
     368             : 
     369             :          ! Evaluate (by interpolation) the norm of the rotational force in the
     370             :          ! minimum of the rotational search (this is for print-out only)
     371        2106 :          ALLOCATE (work(SIZE(dimer_env%nvec)))
     372       24714 :          work = dimer_env%rot%g1
     373             :          work = SIN(dimer_env%rot%angle1 - dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1 + &
     374             :                 SIN(dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1p + &
     375             :                 (1.0_dp - COS(dimer_env%rot%angle2) - SIN(dimer_env%rot%angle2)*TAN(dimer_env%rot%angle1/2.0_dp))* &
     376       24714 :                 dimer_env%rot%g0
     377       24714 :          work = -2.0_dp*(work - dimer_env%rot%g0)
     378       36720 :          work = work - DOT_PRODUCT(work, dimer_env%nvec)*dimer_env%nvec
     379       12708 :          opt_energy = SQRT(DOT_PRODUCT(work, work))
     380         702 :          DEALLOCATE (work)
     381             :       END IF
     382         726 :       dimer_env%rot%angle2 = angle
     383         726 :       CALL timestop(handle)
     384             : 
     385         726 :    END SUBROUTINE rotmin_2pnt
     386             : 
     387             : ! **************************************************************************************************
     388             : !> \brief Line Minimization - Fit
     389             : !> \param gopt_env ...
     390             : !> \param xvec ...
     391             : !> \param xi ...
     392             : !> \param opt_energy ...
     393             : !> \param brack_limit ...
     394             : !> \param step ...
     395             : !> \param output_unit ...
     396             : !> \param gopt_param ...
     397             : !> \param globenv ...
     398             : !> \par History
     399             : !>      10.2005 created [tlaino]
     400             : !> \author Teodoro Laino
     401             : !> \note
     402             : !>      Given as input the vector XVEC and XI, finds the scalar
     403             : !>      xmin that minimizes the energy XVEC+xmin*XI. Replace step
     404             : !>      with the optimal value. Enhanced Version
     405             : ! **************************************************************************************************
     406          54 :    SUBROUTINE linmin_fit(gopt_env, xvec, xi, opt_energy, &
     407             :                          brack_limit, step, output_unit, gopt_param, globenv)
     408             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     409             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xvec, xi
     410             :       REAL(KIND=dp)                                      :: opt_energy, brack_limit, step
     411             :       INTEGER                                            :: output_unit
     412             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     413             :       TYPE(global_environment_type), POINTER             :: globenv
     414             : 
     415             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'linmin_fit'
     416             : 
     417             :       INTEGER                                            :: handle, loc_iter, odim
     418             :       LOGICAL                                            :: should_stop
     419             :       REAL(KIND=dp)                                      :: ax, bx, fprev, rms_dr, rms_force, scale, &
     420             :                                                             xmin, xx
     421          54 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     422          54 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hist
     423             : 
     424          54 :       CALL timeset(routineN, handle)
     425             : 
     426          54 :       NULLIFY (pcom, xicom, hist)
     427          54 :       rms_dr = gopt_param%rms_dr
     428          54 :       rms_force = gopt_param%rms_force
     429         162 :       ALLOCATE (pcom(SIZE(xvec)))
     430         162 :       ALLOCATE (xicom(SIZE(xvec)))
     431             : 
     432       12342 :       pcom = xvec
     433       12342 :       xicom = xi
     434       12342 :       xicom = xicom/SQRT(DOT_PRODUCT(xicom, xicom))
     435          54 :       step = step*0.8_dp ! target a little before the minimum for the first point
     436          54 :       ax = 0.0_dp
     437          54 :       xx = step
     438             :       CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
     439          54 :                      histpoint=hist, globenv=globenv)
     440             :       !
     441          54 :       fprev = 0.0_dp
     442         288 :       opt_energy = MINVAL(hist(:, 2))
     443          54 :       odim = SIZE(hist, 1)
     444          54 :       scale = 0.25_dp
     445          54 :       loc_iter = 0
     446         330 :       DO WHILE (ABS(hist(odim, 3)) > rms_force*scale .OR. ABS(hist(odim, 1) - hist(odim - 1, 1)) > scale*rms_dr)
     447         276 :          CALL external_control(should_stop, "LINFIT", globenv=globenv)
     448         276 :          IF (should_stop) EXIT
     449             :          !
     450         276 :          loc_iter = loc_iter + 1
     451         276 :          fprev = opt_energy
     452         276 :          xmin = FindMin(hist(:, 1), hist(:, 2), hist(:, 3))
     453         276 :          CALL reallocate(hist, 1, odim + 1, 1, 3)
     454         276 :          hist(odim + 1, 1) = xmin
     455         276 :          hist(odim + 1, 3) = cg_deval1d(gopt_env, xmin, pcom, xicom, opt_energy)
     456         276 :          hist(odim + 1, 2) = opt_energy
     457         330 :          odim = SIZE(hist, 1)
     458             :       END DO
     459             :       !
     460        6198 :       xicom = xmin*xicom
     461          54 :       step = xmin
     462       12342 :       xvec = xvec + xicom
     463          54 :       DEALLOCATE (pcom)
     464          54 :       DEALLOCATE (xicom)
     465          54 :       DEALLOCATE (hist)
     466          54 :       IF (output_unit > 0) THEN
     467          27 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     468             :          WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
     469          27 :             "***", "FIT LS  - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
     470          27 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     471             :       END IF
     472          54 :       CALL timestop(handle)
     473             : 
     474          54 :    END SUBROUTINE linmin_fit
     475             : 
     476             : ! **************************************************************************************************
     477             : !> \brief Line Minimization routine - GOLD
     478             : !> \param gopt_env ...
     479             : !> \param xvec ...
     480             : !> \param xi ...
     481             : !> \param opt_energy ...
     482             : !> \param brent_tol ...
     483             : !> \param brent_max_iter ...
     484             : !> \param brack_limit ...
     485             : !> \param step ...
     486             : !> \param output_unit ...
     487             : !> \param globenv ...
     488             : !> \par History
     489             : !>      10.2005 created [tlaino]
     490             : !> \author Teodoro Laino
     491             : !> \note
     492             : !>      Given as input the vector XVEC and XI, finds the scalar
     493             : !>      xmin that minimizes the energy XVEC+xmin*XI. Replaces XMIN
     494             : !>      with the optimal value
     495             : ! **************************************************************************************************
     496         182 :    SUBROUTINE linmin_gold(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
     497             :                           brack_limit, step, output_unit, globenv)
     498             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     499             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xvec, xi
     500             :       REAL(KIND=dp)                                      :: opt_energy, brent_tol
     501             :       INTEGER                                            :: brent_max_iter
     502             :       REAL(KIND=dp)                                      :: brack_limit, step
     503             :       INTEGER                                            :: output_unit
     504             :       TYPE(global_environment_type), POINTER             :: globenv
     505             : 
     506             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'linmin_gold'
     507             : 
     508             :       INTEGER                                            :: handle
     509             :       REAL(KIND=dp)                                      :: ax, bx, xmin, xx
     510         182 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     511             : 
     512         182 :       CALL timeset(routineN, handle)
     513             : 
     514         182 :       NULLIFY (pcom, xicom)
     515         546 :       ALLOCATE (pcom(SIZE(xvec)))
     516         546 :       ALLOCATE (xicom(SIZE(xvec)))
     517             : 
     518       60854 :       pcom = xvec
     519       60854 :       xicom = xi
     520       60854 :       xicom = xicom/SQRT(DOT_PRODUCT(xicom, xicom))
     521         182 :       step = step*0.8_dp ! target a little before the minimum for the first point
     522         182 :       ax = 0.0_dp
     523         182 :       xx = step
     524             :       CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
     525         182 :                      globenv=globenv)
     526             : 
     527             :       opt_energy = cg_dbrent(gopt_env, ax, xx, bx, brent_tol, brent_max_iter, &
     528         182 :                              xmin, pcom, xicom, output_unit, globenv)
     529       30518 :       xicom = xmin*xicom
     530         182 :       step = xmin
     531       60854 :       xvec = xvec + xicom
     532         182 :       DEALLOCATE (pcom)
     533         182 :       DEALLOCATE (xicom)
     534         182 :       CALL timestop(handle)
     535         182 :    END SUBROUTINE linmin_gold
     536             : 
     537             : ! **************************************************************************************************
     538             : !> \brief Routine for initially bracketing a minimum based on the golden search
     539             : !>      minimum
     540             : !> \param gopt_env ...
     541             : !> \param ax ...
     542             : !> \param bx ...
     543             : !> \param cx ...
     544             : !> \param pcom ...
     545             : !> \param xicom ...
     546             : !> \param brack_limit ...
     547             : !> \param output_unit ...
     548             : !> \param histpoint ...
     549             : !> \param globenv ...
     550             : !> \par History
     551             : !>      10.2005 created [tlaino]
     552             : !> \author Teodoro Laino
     553             : !> \note
     554             : !>      Given two distinct initial points ax and bx this routine searches
     555             : !>      in the downhill direction and returns new points ax, bx, cx that
     556             : !>      bracket the minimum of the function
     557             : ! **************************************************************************************************
     558         236 :    SUBROUTINE cg_mnbrak(gopt_env, ax, bx, cx, pcom, xicom, brack_limit, output_unit, &
     559             :                         histpoint, globenv)
     560             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     561             :       REAL(KIND=dp)                                      :: ax, bx, cx
     562             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     563             :       REAL(KIND=dp)                                      :: brack_limit
     564             :       INTEGER                                            :: output_unit
     565             :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: histpoint
     566             :       TYPE(global_environment_type), POINTER             :: globenv
     567             : 
     568             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cg_mnbrak'
     569             : 
     570             :       INTEGER                                            :: handle, loc_iter, odim
     571             :       LOGICAL                                            :: hist, should_stop
     572             :       REAL(KIND=dp)                                      :: dum, fa, fb, fc, fu, gold, q, r, u, ulim
     573             : 
     574         236 :       CALL timeset(routineN, handle)
     575         236 :       hist = PRESENT(histpoint)
     576         236 :       IF (hist) THEN
     577          54 :          CPASSERT(.NOT. ASSOCIATED(histpoint))
     578          54 :          ALLOCATE (histpoint(3, 3))
     579             :       END IF
     580         236 :       gold = (1.0_dp + SQRT(5.0_dp))/2.0_dp
     581         236 :       IF (hist) THEN
     582          54 :          histpoint(1, 1) = ax
     583          54 :          histpoint(1, 3) = cg_deval1d(gopt_env, ax, pcom, xicom, fa)
     584          54 :          histpoint(1, 2) = fa
     585          54 :          histpoint(2, 1) = bx
     586          54 :          histpoint(2, 3) = cg_deval1d(gopt_env, bx, pcom, xicom, fb)
     587          54 :          histpoint(2, 2) = fb
     588             :       ELSE
     589         182 :          fa = cg_eval1d(gopt_env, ax, pcom, xicom)
     590         182 :          fb = cg_eval1d(gopt_env, bx, pcom, xicom)
     591             :       END IF
     592         236 :       IF (fb .GT. fa) THEN
     593          68 :          dum = ax
     594          68 :          ax = bx
     595          68 :          bx = dum
     596          68 :          dum = fb
     597          68 :          fb = fa
     598          68 :          fa = dum
     599             :       END IF
     600         236 :       cx = bx + gold*(bx - ax)
     601         236 :       IF (hist) THEN
     602          54 :          histpoint(3, 1) = cx
     603          54 :          histpoint(3, 3) = cg_deval1d(gopt_env, cx, pcom, xicom, fc)
     604          54 :          histpoint(3, 2) = fc
     605             :       ELSE
     606         182 :          fc = cg_eval1d(gopt_env, cx, pcom, xicom)
     607             :       END IF
     608         236 :       loc_iter = 3
     609         256 :       DO WHILE (fb .GE. fc)
     610          46 :          CALL external_control(should_stop, "MNBRACK", globenv=globenv)
     611          46 :          IF (should_stop) EXIT
     612             :          !
     613          46 :          r = (bx - ax)*(fb - fc)
     614          46 :          q = (bx - cx)*(fb - fa)
     615          46 :          u = bx - ((bx - cx)*q - (bx - ax)*r)/(2.0_dp*SIGN(MAX(ABS(q - r), TINY(0.0_dp)), q - r))
     616          46 :          ulim = bx + brack_limit*(cx - bx)
     617          46 :          IF ((bx - u)*(u - cx) .GT. 0.0_dp) THEN
     618          26 :             IF (hist) THEN
     619          10 :                odim = SIZE(histpoint, 1)
     620          10 :                CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     621          10 :                histpoint(odim + 1, 1) = u
     622          10 :                histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     623          10 :                histpoint(odim + 1, 2) = fu
     624             :             ELSE
     625          16 :                fu = cg_eval1d(gopt_env, u, pcom, xicom)
     626             :             END IF
     627          26 :             loc_iter = loc_iter + 1
     628          26 :             IF (fu .LT. fc) THEN
     629          26 :                ax = bx
     630          26 :                fa = fb
     631          26 :                bx = u
     632          26 :                fb = fu
     633          26 :                EXIT
     634           0 :             ELSE IF (fu .GT. fb) THEN
     635           0 :                cx = u
     636           0 :                fc = fu
     637           0 :                EXIT
     638             :             END IF
     639           0 :             u = cx + gold*(cx - bx)
     640           0 :             IF (hist) THEN
     641           0 :                odim = SIZE(histpoint, 1)
     642           0 :                CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     643           0 :                histpoint(odim + 1, 1) = u
     644           0 :                histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     645           0 :                histpoint(odim + 1, 2) = fu
     646             :             ELSE
     647           0 :                fu = cg_eval1d(gopt_env, u, pcom, xicom)
     648             :             END IF
     649           0 :             loc_iter = loc_iter + 1
     650          20 :          ELSE IF ((cx - u)*(u - ulim) .GT. 0.) THEN
     651          20 :             IF (hist) THEN
     652           4 :                odim = SIZE(histpoint, 1)
     653           4 :                CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     654           4 :                histpoint(odim + 1, 1) = u
     655           4 :                histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     656           4 :                histpoint(odim + 1, 2) = fu
     657             :             ELSE
     658          16 :                fu = cg_eval1d(gopt_env, u, pcom, xicom)
     659             :             END IF
     660          20 :             loc_iter = loc_iter + 1
     661          20 :             IF (fu .LT. fc) THEN
     662          18 :                bx = cx
     663          18 :                cx = u
     664          18 :                u = cx + gold*(cx - bx)
     665          18 :                fb = fc
     666          18 :                fc = fu
     667          18 :                IF (hist) THEN
     668           4 :                   odim = SIZE(histpoint, 1)
     669           4 :                   CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     670           4 :                   histpoint(odim + 1, 1) = u
     671           4 :                   histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     672           4 :                   histpoint(odim + 1, 2) = fu
     673             :                ELSE
     674          14 :                   fu = cg_eval1d(gopt_env, u, pcom, xicom)
     675             :                END IF
     676          18 :                loc_iter = loc_iter + 1
     677             :             END IF
     678           0 :          ELSE IF ((u - ulim)*(ulim - cx) .GE. 0.) THEN
     679           0 :             u = ulim
     680           0 :             IF (hist) THEN
     681           0 :                odim = SIZE(histpoint, 1)
     682           0 :                CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     683           0 :                histpoint(odim + 1, 1) = u
     684           0 :                histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     685           0 :                histpoint(odim + 1, 2) = fu
     686             :             ELSE
     687           0 :                fu = cg_eval1d(gopt_env, u, pcom, xicom)
     688             :             END IF
     689           0 :             loc_iter = loc_iter + 1
     690             :          ELSE
     691           0 :             u = cx + gold*(cx - bx)
     692           0 :             IF (hist) THEN
     693           0 :                odim = SIZE(histpoint, 1)
     694           0 :                CALL reallocate(histpoint, 1, odim + 1, 1, 3)
     695           0 :                histpoint(odim + 1, 1) = u
     696           0 :                histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     697           0 :                histpoint(odim + 1, 2) = fu
     698             :             ELSE
     699           0 :                fu = cg_eval1d(gopt_env, u, pcom, xicom)
     700             :             END IF
     701           0 :             loc_iter = loc_iter + 1
     702             :          END IF
     703          20 :          ax = bx
     704          20 :          bx = cx
     705          20 :          cx = u
     706          20 :          fa = fb
     707          20 :          fb = fc
     708          20 :          fc = fu
     709             :       END DO
     710         236 :       IF (output_unit > 0) THEN
     711         118 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     712             :          WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
     713         118 :             "***", "MNBRACK - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
     714         118 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     715             :       END IF
     716         236 :       CALL timestop(handle)
     717         236 :    END SUBROUTINE cg_mnbrak
     718             : 
     719             : ! **************************************************************************************************
     720             : !> \brief Routine implementing the Brent Method
     721             : !>      Brent,R.P. Algorithm for Minimization without Derivatives, Chapt.5
     722             : !>      1973
     723             : !>      Extension in the use of derivatives
     724             : !> \param gopt_env ...
     725             : !> \param ax ...
     726             : !> \param bx ...
     727             : !> \param cx ...
     728             : !> \param tol ...
     729             : !> \param itmax ...
     730             : !> \param xmin ...
     731             : !> \param pcom ...
     732             : !> \param xicom ...
     733             : !> \param output_unit ...
     734             : !> \param globenv ...
     735             : !> \return ...
     736             : !> \par History
     737             : !>      10.2005 created [tlaino]
     738             : !> \author Teodoro Laino
     739             : !> \note
     740             : !>      Given a bracketing  triplet of abscissas ax, bx, cx (such that bx
     741             : !>      is between ax and cx and energy of bx is less than energy of ax and cx),
     742             : !>      this routine isolates the minimum to a precision of about tol using
     743             : !>      Brent method. This routine implements the extension of the Brent Method
     744             : !>      using derivatives
     745             : ! **************************************************************************************************
     746         182 :    FUNCTION cg_dbrent(gopt_env, ax, bx, cx, tol, itmax, xmin, pcom, xicom, output_unit, &
     747             :                       globenv) RESULT(dbrent)
     748             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     749             :       REAL(KIND=dp)                                      :: ax, bx, cx, tol
     750             :       INTEGER                                            :: itmax
     751             :       REAL(KIND=dp)                                      :: xmin
     752             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     753             :       INTEGER                                            :: output_unit
     754             :       TYPE(global_environment_type), POINTER             :: globenv
     755             :       REAL(KIND=dp)                                      :: dbrent
     756             : 
     757             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cg_dbrent'
     758             :       REAL(KIND=dp), PARAMETER                           :: zeps = 1.0E-8_dp
     759             : 
     760             :       INTEGER                                            :: handle, iter, loc_iter
     761             :       LOGICAL                                            :: ok1, ok2, should_stop, skip0, skip1
     762             :       REAL(KIND=dp)                                      :: a, b, d, d1, d2, du, dv, dw, dx, e, fu, &
     763             :                                                             fv, fw, fx, olde, tol1, tol2, u, u1, &
     764             :                                                             u2, v, w, x, xm
     765             : 
     766         182 :       CALL timeset(routineN, handle)
     767         182 :       a = MIN(ax, cx)
     768         182 :       b = MAX(ax, cx)
     769         182 :       v = bx; w = v; x = v
     770         182 :       e = 0.0_dp
     771         182 :       dx = cg_deval1d(gopt_env, x, pcom, xicom, fx)
     772         182 :       fv = fx
     773         182 :       fw = fx
     774         182 :       dv = dx
     775         182 :       dw = dx
     776         182 :       loc_iter = 1
     777         646 :       DO iter = 1, itmax
     778         646 :          CALL external_control(should_stop, "BRENT", globenv=globenv)
     779         646 :          IF (should_stop) EXIT
     780             :          !
     781         646 :          xm = 0.5_dp*(a + b)
     782         646 :          tol1 = tol*ABS(x) + zeps
     783         646 :          tol2 = 2.0_dp*tol1
     784         646 :          skip0 = .FALSE.
     785         646 :          skip1 = .FALSE.
     786         646 :          IF (ABS(x - xm) .LE. (tol2 - 0.5_dp*(b - a))) EXIT
     787         610 :          IF (ABS(e) .GT. tol1) THEN
     788         410 :             d1 = 2.0_dp*(b - a)
     789         410 :             d2 = d1
     790         410 :             IF (dw .NE. dx) d1 = (w - x)*dx/(dx - dw)
     791         410 :             IF (dv .NE. dx) d2 = (v - x)*dx/(dx - dv)
     792         410 :             u1 = x + d1
     793         410 :             u2 = x + d2
     794         410 :             ok1 = ((a - u1)*(u1 - b) .GT. 0.0_dp) .AND. (dx*d1 .LE. 0.0_dp)
     795         410 :             ok2 = ((a - u2)*(u2 - b) .GT. 0.0_dp) .AND. (dx*d2 .LE. 0.0_dp)
     796         410 :             olde = e
     797         410 :             e = d
     798         410 :             IF (.NOT. (ok1 .OR. ok2)) THEN
     799             :                skip0 = .TRUE.
     800         394 :             ELSE IF (ok1 .AND. ok2) THEN
     801         294 :                IF (ABS(d1) .LT. ABS(d2)) THEN
     802             :                   d = d1
     803             :                ELSE
     804         172 :                   d = d2
     805             :                END IF
     806         100 :             ELSE IF (ok1) THEN
     807             :                d = d1
     808             :             ELSE
     809           0 :                d = d2
     810             :             END IF
     811         410 :             IF (.NOT. skip0) THEN
     812         394 :                IF (ABS(d) .GT. ABS(0.5_dp*olde)) skip0 = .TRUE.
     813             :                IF (.NOT. skip0) THEN
     814         392 :                   u = x + d
     815         392 :                   IF ((u - a) .LT. tol2 .OR. (b - u) .LT. tol2) d = SIGN(tol1, xm - x)
     816             :                   skip1 = .TRUE.
     817             :                END IF
     818             :             END IF
     819             :          END IF
     820         610 :          IF (.NOT. skip1) THEN
     821         218 :             IF (dx .GE. 0.0_dp) THEN
     822         104 :                e = a - x
     823             :             ELSE
     824         114 :                e = b - x
     825             :             END IF
     826         218 :             d = 0.5_dp*e
     827             :          END IF
     828         610 :          IF (ABS(d) .GE. tol1) THEN
     829         426 :             u = x + d
     830         426 :             du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     831         426 :             loc_iter = loc_iter + 1
     832             :          ELSE
     833         184 :             u = x + SIGN(tol1, d)
     834         184 :             du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
     835         184 :             loc_iter = loc_iter + 1
     836         184 :             IF (fu .GT. fx) EXIT
     837             :          END IF
     838        1292 :          IF (fu .LE. fx) THEN
     839         316 :             IF (u .GE. x) THEN
     840             :                a = x
     841             :             ELSE
     842         132 :                b = x
     843             :             END IF
     844         316 :             v = w; fv = fw; dv = dw; w = x
     845         316 :             fw = fx; dw = dx; x = u; fx = fu; dx = du
     846             :          ELSE
     847         148 :             IF (u .LT. x) THEN
     848             :                a = u
     849             :             ELSE
     850          76 :                b = u
     851             :             END IF
     852         148 :             IF (fu .LE. fw .OR. w .EQ. x) THEN
     853             :                v = w; fv = fw; dv = dw
     854             :                w = u; fw = fu; dw = du
     855          14 :             ELSE IF (fu .LE. fv .OR. v .EQ. x .OR. v .EQ. w) THEN
     856          14 :                v = u
     857          14 :                fv = fu
     858          14 :                dv = du
     859             :             END IF
     860             :          END IF
     861             :       END DO
     862         182 :       IF (output_unit > 0) THEN
     863          91 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     864             :          WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
     865          91 :             "***", "BRENT   - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
     866          91 :          IF (iter == itmax + 1) &
     867             :             WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,T78,A)") &
     868           0 :             "***", "BRENT - NUMBER OF ITERATIONS EXCEEDED ", "***"
     869          91 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     870             :       END IF
     871         182 :       CPASSERT(iter /= itmax + 1)
     872         182 :       xmin = x
     873         182 :       dbrent = fx
     874         182 :       CALL timestop(handle)
     875             : 
     876         182 :    END FUNCTION cg_dbrent
     877             : 
     878             : ! **************************************************************************************************
     879             : !> \brief Evaluates energy in one dimensional space defined by the point
     880             : !>      pcom and with direction xicom, position x
     881             : !> \param gopt_env ...
     882             : !> \param x ...
     883             : !> \param pcom ...
     884             : !> \param xicom ...
     885             : !> \return ...
     886             : !> \par History
     887             : !>      10.2005 created [tlaino]
     888             : !> \author Teodoro Laino
     889             : ! **************************************************************************************************
     890         592 :    FUNCTION cg_eval1d(gopt_env, x, pcom, xicom) RESULT(my_val)
     891             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     892             :       REAL(KIND=dp)                                      :: x
     893             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     894             :       REAL(KIND=dp)                                      :: my_val
     895             : 
     896             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cg_eval1d'
     897             : 
     898             :       INTEGER                                            :: handle
     899         592 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: xvec
     900             : 
     901         592 :       CALL timeset(routineN, handle)
     902             : 
     903        1776 :       ALLOCATE (xvec(SIZE(pcom)))
     904      205868 :       xvec = pcom + x*xicom
     905             :       CALL cp_eval_at(gopt_env, xvec, my_val, master=gopt_env%force_env%para_env%mepos, &
     906         592 :                       final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     907         592 :       DEALLOCATE (xvec)
     908             : 
     909         592 :       CALL timestop(handle)
     910             : 
     911         592 :    END FUNCTION cg_eval1d
     912             : 
     913             : ! **************************************************************************************************
     914             : !> \brief Evaluates derivatives in one dimensional space defined by the point
     915             : !>      pcom and with direction xicom, position x
     916             : !> \param gopt_env ...
     917             : !> \param x ...
     918             : !> \param pcom ...
     919             : !> \param xicom ...
     920             : !> \param fval ...
     921             : !> \return ...
     922             : !> \par History
     923             : !>      10.2005 created [tlaino]
     924             : !> \author Teodoro Laino
     925             : ! **************************************************************************************************
     926        1248 :    FUNCTION cg_deval1d(gopt_env, x, pcom, xicom, fval) RESULT(my_val)
     927             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     928             :       REAL(KIND=dp)                                      :: x
     929             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pcom, xicom
     930             :       REAL(KIND=dp)                                      :: fval, my_val
     931             : 
     932             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cg_deval1d'
     933             : 
     934             :       INTEGER                                            :: handle
     935             :       REAL(KIND=dp)                                      :: energy
     936        1248 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: grad, xvec
     937             : 
     938        1248 :       CALL timeset(routineN, handle)
     939             : 
     940        3744 :       ALLOCATE (xvec(SIZE(pcom)))
     941        3744 :       ALLOCATE (grad(SIZE(pcom)))
     942      354168 :       xvec = pcom + x*xicom
     943             :       CALL cp_eval_at(gopt_env, xvec, energy, grad, master=gopt_env%force_env%para_env%mepos, &
     944        1248 :                       final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
     945      177084 :       my_val = DOT_PRODUCT(grad, xicom)
     946        1248 :       fval = energy
     947        1248 :       DEALLOCATE (xvec)
     948        1248 :       DEALLOCATE (grad)
     949        1248 :       CALL timestop(handle)
     950             : 
     951        1248 :    END FUNCTION cg_deval1d
     952             : 
     953             : ! **************************************************************************************************
     954             : !> \brief Find the minimum of a parabolic function obtained with a least square fit
     955             : !> \param x ...
     956             : !> \param y ...
     957             : !> \param dy ...
     958             : !> \return ...
     959             : !> \par History
     960             : !>      10.2005 created [fawzi]
     961             : !> \author Fawzi Mohamed
     962             : ! **************************************************************************************************
     963         276 :    FUNCTION FindMin(x, y, dy) RESULT(res)
     964             :       REAL(kind=dp), DIMENSION(:)                        :: x, y, dy
     965             :       REAL(kind=dp)                                      :: res
     966             : 
     967             :       INTEGER                                            :: i, info, iwork(8*3), lwork, min_pos, np
     968             :       REAL(kind=dp)                                      :: diag(3), res1(3), res2(3), res3(3), &
     969             :                                                             spread, sum_x, sum_xx, tmpw(1), &
     970             :                                                             vt(3, 3)
     971         276 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: work
     972         552 :       REAL(kind=dp), DIMENSION(2*SIZE(x), 3)             :: f
     973         552 :       REAL(kind=dp), DIMENSION(2*SIZE(x))                :: b, w
     974         552 :       REAL(kind=dp)                                      :: u(2*SIZE(x), 3)
     975             : 
     976         276 :       np = SIZE(x)
     977         276 :       CPASSERT(np > 1)
     978         276 :       sum_x = 0._dp
     979         276 :       sum_xx = 0._dp
     980         276 :       min_pos = 1
     981        2178 :       DO i = 1, np
     982        1902 :          sum_xx = sum_xx + x(i)**2
     983        1902 :          sum_x = sum_x + x(i)
     984        2178 :          IF (y(min_pos) > y(i)) min_pos = i
     985             :       END DO
     986             :       spread = SQRT(sum_xx/REAL(np, dp) - (sum_x/REAL(np, dp))**2)
     987        2178 :       DO i = 1, np
     988        1902 :          w(i) = EXP(-(REAL(np - i, dp))**2/(REAL(2*9, dp)))
     989        2178 :          w(i + np) = 2._dp*w(i)
     990             :       END DO
     991        2178 :       DO i = 1, np
     992        1902 :          f(i, 1) = w(i)
     993        1902 :          f(i, 2) = x(i)*w(i)
     994        1902 :          f(i, 3) = x(i)**2*w(i)
     995        1902 :          f(i + np, 1) = 0
     996        1902 :          f(i + np, 2) = w(i + np)
     997        2178 :          f(i + np, 3) = 2*x(i)*w(i + np)
     998             :       END DO
     999        2178 :       DO i = 1, np
    1000        1902 :          b(i) = y(i)*w(i)
    1001        2178 :          b(i + np) = dy(i)*w(i + np)
    1002             :       END DO
    1003         276 :       lwork = -1
    1004             :       CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), tmpw, lwork, &
    1005         276 :                   iwork, info)
    1006         276 :       lwork = CEILING(tmpw(1))
    1007         828 :       ALLOCATE (work(lwork))
    1008             :       CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), work, lwork, &
    1009         276 :                   iwork, info)
    1010         276 :       DEALLOCATE (work)
    1011         276 :       CALL dgemv('T', SIZE(u, 1), SIZE(u, 2), 1._dp, u, SIZE(u, 1), b, 1, 0._dp, res1, 1)
    1012        1104 :       DO i = 1, 3
    1013        1104 :          res2(i) = res1(i)/diag(i)
    1014             :       END DO
    1015         276 :       CALL dgemv('T', 3, 3, 1._dp, vt, SIZE(vt, 1), res2, 1, 0._dp, res3, 1)
    1016         276 :       res = -0.5*res3(2)/res3(3)
    1017         276 :    END FUNCTION FindMin
    1018             : 
    1019             : ! **************************************************************************************************
    1020             : !> \brief Computes the Conjugate direction for the next search
    1021             : !> \param gopt_env ...
    1022             : !> \param Fletcher_Reeves ...
    1023             : !> \param g contains the theta  of the previous step.. (norm 1.0 vector)
    1024             : !> \param xi contains the -theta of the present step.. (norm 1.0 vector)
    1025             : !> \param h contains the search direction of the previous step (must be orthogonal
    1026             : !>            to nvec of the previous step (nvec_old))
    1027             : !> \par   Info for DIMER method
    1028             : !> \par History
    1029             : !>      10.2005 created [tlaino]
    1030             : !> \author Teodoro Laino
    1031             : ! **************************************************************************************************
    1032        2052 :    SUBROUTINE get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
    1033             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
    1034             :       LOGICAL, INTENT(IN)                                :: Fletcher_Reeves
    1035             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: g, xi, h
    1036             : 
    1037             :       CHARACTER(len=*), PARAMETER :: routineN = 'get_conjugate_direction'
    1038             : 
    1039             :       INTEGER                                            :: handle
    1040             :       LOGICAL                                            :: check
    1041             :       REAL(KIND=dp)                                      :: dgg, gam, gg, norm, norm_h
    1042             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
    1043             : 
    1044        2052 :       CALL timeset(routineN, handle)
    1045        2052 :       NULLIFY (dimer_env)
    1046        2052 :       IF (.NOT. gopt_env%dimer_rotation) THEN
    1047      311820 :          gg = DOT_PRODUCT(g, g)
    1048        1458 :          IF (Fletcher_Reeves) THEN
    1049           0 :             dgg = DOT_PRODUCT(xi, xi)
    1050             :          ELSE
    1051      311820 :             dgg = DOT_PRODUCT((xi + g), xi)
    1052             :          END IF
    1053        1458 :          gam = dgg/gg
    1054      622182 :          g = h
    1055      622182 :          h = -xi + gam*h
    1056             :       ELSE
    1057         594 :          dimer_env => gopt_env%dimer_env
    1058       10932 :          check = ABS(DOT_PRODUCT(g, g) - 1.0_dp) < MAX(1.0E-9_dp, dimer_thrs)
    1059         594 :          CPASSERT(check)
    1060             : 
    1061       10932 :          check = ABS(DOT_PRODUCT(xi, xi) - 1.0_dp) < MAX(1.0E-9_dp, dimer_thrs)
    1062         594 :          CPASSERT(check)
    1063             : 
    1064       10932 :          check = ABS(DOT_PRODUCT(h, dimer_env%cg_rot%nvec_old)) < MAX(1.0E-9_dp, dimer_thrs)
    1065         594 :          CPASSERT(check)
    1066         594 :          gg = dimer_env%cg_rot%norm_theta_old**2
    1067         594 :          IF (Fletcher_Reeves) THEN
    1068           0 :             dgg = dimer_env%cg_rot%norm_theta**2
    1069             :          ELSE
    1070         594 :             norm = dimer_env%cg_rot%norm_theta*dimer_env%cg_rot%norm_theta_old
    1071       10932 :             dgg = dimer_env%cg_rot%norm_theta**2 + DOT_PRODUCT(g, xi)*norm
    1072             :          END IF
    1073             :          ! Compute Theta** and store it in nvec_old
    1074         594 :          CALL rotate_dimer(dimer_env%cg_rot%nvec_old, g, dimer_env%rot%angle2 + pi/2.0_dp)
    1075         594 :          gam = dgg/gg
    1076       21270 :          g = h
    1077       21270 :          h = -xi*dimer_env%cg_rot%norm_theta + gam*dimer_env%cg_rot%norm_h*dimer_env%cg_rot%nvec_old
    1078       31608 :          h = h - DOT_PRODUCT(h, dimer_env%nvec)*dimer_env%nvec
    1079       10932 :          norm_h = SQRT(DOT_PRODUCT(h, h))
    1080         594 :          IF (norm_h < EPSILON(0.0_dp)) THEN
    1081           0 :             h = 0.0_dp
    1082             :          ELSE
    1083       10932 :             h = h/norm_h
    1084             :          END IF
    1085         594 :          dimer_env%cg_rot%norm_h = norm_h
    1086             :       END IF
    1087        2052 :       CALL timestop(handle)
    1088             : 
    1089        2052 :    END SUBROUTINE get_conjugate_direction
    1090             : 
    1091             : END MODULE cg_utils

Generated by: LCOV version 1.15