Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Optimizers used by pao_main.F
      10             : !> \author Ole Schuett
      11             : ! **************************************************************************************************
      12             : MODULE pao_optimizer
      13             :    USE arnoldi_api,                     ONLY: arnoldi_extremal
      14             :    USE cp_dbcsr_api,                    ONLY: &
      15             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_release, &
      16             :         dbcsr_scale, dbcsr_set, dbcsr_type
      17             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_add_on_diag,&
      18             :                                               dbcsr_dot,&
      19             :                                               dbcsr_frobenius_norm,&
      20             :                                               dbcsr_reserve_diag_blocks
      21             :    USE kinds,                           ONLY: dp
      22             :    USE pao_input,                       ONLY: pao_opt_bfgs,&
      23             :                                               pao_opt_cg
      24             :    USE pao_types,                       ONLY: pao_env_type
      25             : #include "./base/base_uses.f90"
      26             : 
      27             :    IMPLICIT NONE
      28             : 
      29             :    PRIVATE
      30             : 
      31             :    PUBLIC :: pao_opt_init, pao_opt_finalize, pao_opt_new_dir
      32             : 
      33             : CONTAINS
      34             : 
      35             : ! **************************************************************************************************
      36             : !> \brief Initialize the optimizer
      37             : !> \param pao ...
      38             : ! **************************************************************************************************
      39         246 :    SUBROUTINE pao_opt_init(pao)
      40             :       TYPE(pao_env_type), POINTER                        :: pao
      41             : 
      42         246 :       CALL dbcsr_copy(pao%matrix_D, pao%matrix_G)
      43         246 :       CALL dbcsr_set(pao%matrix_D, 0.0_dp)
      44             : 
      45         246 :       CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_D)
      46             : 
      47         246 :       IF (pao%precondition) THEN
      48          82 :          CALL dbcsr_copy(pao%matrix_D_preconed, pao%matrix_D)
      49             :       END IF
      50             : 
      51         246 :       IF (pao%optimizer == pao_opt_bfgs) &
      52          12 :          CALL pao_opt_init_bfgs(pao)
      53             : 
      54         246 :    END SUBROUTINE pao_opt_init
      55             : 
      56             : ! **************************************************************************************************
      57             : !> \brief Initialize the BFGS optimizer
      58             : !> \param pao ...
      59             : ! **************************************************************************************************
      60          12 :    SUBROUTINE pao_opt_init_bfgs(pao)
      61             :       TYPE(pao_env_type), POINTER                        :: pao
      62             : 
      63          12 :       INTEGER, DIMENSION(:), POINTER                     :: nparams
      64             : 
      65          12 :       CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nparams)
      66             : 
      67             :       CALL dbcsr_create(pao%matrix_BFGS, &
      68             :                         template=pao%matrix_X, &
      69             :                         row_blk_size=nparams, &
      70             :                         col_blk_size=nparams, &
      71          12 :                         name="PAO matrix_BFGS")
      72             : 
      73          12 :       CALL dbcsr_reserve_diag_blocks(pao%matrix_BFGS)
      74          12 :       CALL dbcsr_set(pao%matrix_BFGS, 0.0_dp)
      75          12 :       CALL dbcsr_add_on_diag(pao%matrix_BFGS, 1.0_dp)
      76             : 
      77          12 :    END SUBROUTINE pao_opt_init_bfgs
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief Finalize the optimizer
      81             : !> \param pao ...
      82             : ! **************************************************************************************************
      83         246 :    SUBROUTINE pao_opt_finalize(pao)
      84             :       TYPE(pao_env_type), POINTER                        :: pao
      85             : 
      86         246 :       CALL dbcsr_release(pao%matrix_D)
      87         246 :       CALL dbcsr_release(pao%matrix_G_prev)
      88         246 :       IF (pao%precondition) &
      89          82 :          CALL dbcsr_release(pao%matrix_D_preconed)
      90             : 
      91         246 :       IF (pao%optimizer == pao_opt_bfgs) &
      92          12 :          CALL dbcsr_release(pao%matrix_BFGS)
      93             : 
      94         246 :    END SUBROUTINE pao_opt_finalize
      95             : 
      96             : ! **************************************************************************************************
      97             : !> \brief Calculates the new search direction.
      98             : !> \param pao ...
      99             : !> \param icycle ...
     100             : ! **************************************************************************************************
     101        2616 :    SUBROUTINE pao_opt_new_dir(pao, icycle)
     102             :       TYPE(pao_env_type), POINTER                        :: pao
     103             :       INTEGER, INTENT(IN)                                :: icycle
     104             : 
     105             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pao_opt_new_dir'
     106             : 
     107             :       INTEGER                                            :: handle
     108             :       TYPE(dbcsr_type)                                   :: matrix_G_preconed
     109             : 
     110        2616 :       CALL timeset(routineN, handle)
     111             : 
     112        2616 :       IF (pao%precondition) THEN
     113             :          ! We can't convert matrix_D for and back every time, the numeric noise would disturb the CG,
     114             :          ! hence we keep matrix_D_preconed around.
     115        1290 :          CALL dbcsr_copy(matrix_G_preconed, pao%matrix_G)
     116             :          CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_G, &
     117        1290 :                              0.0_dp, matrix_G_preconed, retain_sparsity=.TRUE.)
     118        1290 :          CALL pao_opt_new_dir_low(pao, icycle, matrix_G_preconed, pao%matrix_G_prev, pao%matrix_D_preconed)
     119             :          CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_D_preconed, &
     120        1290 :                              0.0_dp, pao%matrix_D, retain_sparsity=.TRUE.)
     121             : 
     122             :          ! store preconditioned gradient for next iteration
     123        1290 :          CALL dbcsr_copy(pao%matrix_G_prev, matrix_G_preconed)
     124             : 
     125        1290 :          pao%norm_G = dbcsr_frobenius_norm(matrix_G_preconed)
     126        1290 :          IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of preconditioned gradient:", pao%norm_G
     127        1290 :          CALL dbcsr_release(matrix_G_preconed)
     128             : 
     129             :       ELSE
     130        1326 :          CALL pao_opt_new_dir_low(pao, icycle, pao%matrix_G, pao%matrix_G_prev, pao%matrix_D)
     131        1326 :          CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_G) ! store gradient for next iteration
     132        1326 :          pao%norm_G = dbcsr_frobenius_norm(pao%matrix_G)
     133        1326 :          IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of gradient:", pao%norm_G
     134             :       END IF
     135             : 
     136        2616 :       CALL timestop(handle)
     137             : 
     138        2616 :    END SUBROUTINE pao_opt_new_dir
     139             : 
     140             : ! **************************************************************************************************
     141             : !> \brief Calculates the new search direction.
     142             : !> \param pao ...
     143             : !> \param icycle ...
     144             : !> \param matrix_G ...
     145             : !> \param matrix_G_prev ...
     146             : !> \param matrix_D ...
     147             : ! **************************************************************************************************
     148        2616 :    SUBROUTINE pao_opt_new_dir_low(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
     149             :       TYPE(pao_env_type), POINTER                        :: pao
     150             :       INTEGER, INTENT(IN)                                :: icycle
     151             :       TYPE(dbcsr_type)                                   :: matrix_G, matrix_G_prev, matrix_D
     152             : 
     153        4944 :       SELECT CASE (pao%optimizer)
     154             :       CASE (pao_opt_cg)
     155        2328 :          CALL pao_opt_newdir_cg(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
     156             :       CASE (pao_opt_bfgs)
     157         288 :          CALL pao_opt_newdir_bfgs(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
     158             :       CASE DEFAULT
     159        2616 :          CPABORT("PAO: unknown optimizer")
     160             :       END SELECT
     161             : 
     162        2616 :    END SUBROUTINE pao_opt_new_dir_low
     163             : 
     164             : ! **************************************************************************************************
     165             : !> \brief Conjugate Gradient algorithm
     166             : !> \param pao ...
     167             : !> \param icycle ...
     168             : !> \param matrix_G ...
     169             : !> \param matrix_G_prev ...
     170             : !> \param matrix_D ...
     171             : ! **************************************************************************************************
     172        2328 :    SUBROUTINE pao_opt_newdir_cg(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
     173             :       TYPE(pao_env_type), POINTER                        :: pao
     174             :       INTEGER, INTENT(IN)                                :: icycle
     175             :       TYPE(dbcsr_type)                                   :: matrix_G, matrix_G_prev, matrix_D
     176             : 
     177             :       REAL(KIND=dp)                                      :: beta, change, trace_D, trace_D_Gnew, &
     178             :                                                             trace_G_mix, trace_G_new, trace_G_prev
     179             : 
     180             :       ! determine CG mixing factor
     181        2328 :       IF (icycle <= pao%cg_init_steps) THEN
     182         444 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| warming up with steepest descent"
     183         444 :          beta = 0.0_dp
     184             :       ELSE
     185        1884 :          CALL dbcsr_dot(matrix_G, matrix_G, trace_G_new)
     186        1884 :          CALL dbcsr_dot(matrix_G_prev, matrix_G_prev, trace_G_prev)
     187        1884 :          CALL dbcsr_dot(matrix_G, matrix_G_prev, trace_G_mix)
     188        1884 :          CALL dbcsr_dot(matrix_D, matrix_G, trace_D_Gnew)
     189        1884 :          CALL dbcsr_dot(matrix_D, matrix_D, trace_D)
     190        1884 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_new ", trace_G_new
     191        1884 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_prev ", trace_G_prev
     192        1884 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_mix ", trace_G_mix
     193        1884 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D ", trace_D
     194        1884 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D_Gnew", trace_D_Gnew
     195             : 
     196        1884 :          IF (trace_G_prev /= 0.0_dp) THEN
     197        1884 :             beta = (trace_G_new - trace_G_mix)/trace_G_prev !Polak-Ribiere
     198             :          END IF
     199             : 
     200        1884 :          IF (beta < 0.0_dp) THEN
     201          78 :             IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because beta < 0"
     202          78 :             beta = 0.0_dp
     203             :          END IF
     204             : 
     205        1884 :          change = trace_D_Gnew**2/trace_D*trace_G_new
     206        1884 :          IF (change > pao%cg_reset_limit) THEN
     207           0 :             IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because change > CG_RESET_LIMIT"
     208           0 :             beta = 0.0_dp
     209             :          END IF
     210             : 
     211             :       END IF
     212             : 
     213        2328 :       IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| beta: ", beta
     214             : 
     215             :       ! calculate new CG direction matrix_D
     216        2328 :       CALL dbcsr_add(matrix_D, matrix_G, beta, -1.0_dp)
     217             : 
     218        2328 :    END SUBROUTINE pao_opt_newdir_cg
     219             : 
     220             : ! **************************************************************************************************
     221             : !> \brief Broyden-Fletcher-Goldfarb-Shanno algorithm
     222             : !> \param pao ...
     223             : !> \param icycle ...
     224             : !> \param matrix_G ...
     225             : !> \param matrix_G_prev ...
     226             : !> \param matrix_D ...
     227             : ! **************************************************************************************************
     228         288 :    SUBROUTINE pao_opt_newdir_bfgs(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
     229             :       TYPE(pao_env_type), POINTER                        :: pao
     230             :       INTEGER, INTENT(IN)                                :: icycle
     231             :       TYPE(dbcsr_type)                                   :: matrix_G, matrix_G_prev, matrix_D
     232             : 
     233             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_opt_newdir_bfgs'
     234             : 
     235             :       INTEGER                                            :: handle
     236             :       LOGICAL                                            :: arnoldi_converged
     237             :       REAL(dp)                                           :: eval_max, eval_min, theta, trace_ry, &
     238             :                                                             trace_sy, trace_yHy, trace_yy
     239             :       TYPE(dbcsr_type)                                   :: matrix_Hy, matrix_Hyr, matrix_r, &
     240             :                                                             matrix_rr, matrix_ryH, matrix_ryHyr, &
     241             :                                                             matrix_s, matrix_y, matrix_yr
     242             : 
     243         288 :       CALL timeset(routineN, handle)
     244             : 
     245             :       !TODO add filtering?
     246             : 
     247             :       ! Notation according to the book from Nocedal and Wright, see chapter 6.
     248         288 :       IF (icycle > 1) THEN
     249             :          ! y = G - G_prev
     250         276 :          CALL dbcsr_copy(matrix_y, matrix_G)
     251         276 :          CALL dbcsr_add(matrix_y, matrix_G_prev, 1.0_dp, -1.0_dp) ! dG
     252             : 
     253             :          ! s = X - X_prev
     254         276 :          CALL dbcsr_copy(matrix_s, matrix_D)
     255         276 :          CALL dbcsr_scale(matrix_s, pao%linesearch%step_size) ! dX
     256             : 
     257             :          ! sy = MATMUL(TRANPOSE(s), y)
     258         276 :          CALL dbcsr_dot(matrix_s, matrix_y, trace_sy)
     259             : 
     260             :          ! heuristic initialization
     261         276 :          IF (icycle == 2) THEN
     262          10 :             CALL dbcsr_dot(matrix_Y, matrix_Y, trace_yy)
     263          10 :             CALL dbcsr_scale(pao%matrix_BFGS, trace_sy/trace_yy)
     264          10 :             IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Initializing with:", trace_sy/trace_yy
     265             :          END IF
     266             : 
     267             :          ! Hy = MATMUL(H, y)
     268         276 :          CALL dbcsr_create(matrix_Hy, template=matrix_G, matrix_type="N")
     269         276 :          CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_y, 0.0_dp, matrix_Hy)
     270             : 
     271             :          ! yHy = MATMUL(TRANPOSE(y), Hy)
     272         276 :          CALL dbcsr_dot(matrix_y, matrix_Hy, trace_yHy)
     273             : 
     274             :          ! Use damped BFGS algorithm to ensure H remains positive definite.
     275             :          ! See chapter 18 in Nocedal and Wright's book for details.
     276             :          ! The formulas were adopted to inverse Hessian algorithm.
     277         276 :          IF (trace_sy < 0.2_dp*trace_yHy) THEN
     278           0 :             theta = 0.8_dp*trace_yHy/(trace_yHy - trace_sy)
     279           0 :             IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Dampening theta:", theta
     280             :          ELSE
     281         276 :             theta = 1.0
     282             :          END IF
     283             : 
     284             :          ! r = theta*s + (1-theta)*Hy
     285         276 :          CALL dbcsr_copy(matrix_r, matrix_s)
     286         276 :          CALL dbcsr_add(matrix_r, matrix_Hy, theta, (1.0_dp - theta))
     287             : 
     288             :          ! use t instead of y to update B matrix
     289         276 :          CALL dbcsr_dot(matrix_r, matrix_y, trace_ry)
     290         276 :          CPASSERT(trace_RY > 0.0_dp)
     291             : 
     292             :          ! yr = MATMUL(y, TRANSPOSE(r))
     293         276 :          CALL dbcsr_create(matrix_yr, template=pao%matrix_BFGS, matrix_type="N")
     294         276 :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_y, matrix_r, 0.0_dp, matrix_yr)
     295             : 
     296             :          ! Hyr = MATMUL(H, yr)
     297         276 :          CALL dbcsr_create(matrix_Hyr, template=pao%matrix_BFGS, matrix_type="N")
     298         276 :          CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_yr, 0.0_dp, matrix_Hyr)
     299             : 
     300             :          ! ryH = MATMUL(TRANSPOSE(yr), H)
     301         276 :          CALL dbcsr_create(matrix_ryH, template=pao%matrix_BFGS, matrix_type="N")
     302         276 :          CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_yr, pao%matrix_BFGS, 0.0_dp, matrix_ryH)
     303             : 
     304             :          ! ryHry = MATMUL(ryH,yr)
     305         276 :          CALL dbcsr_create(matrix_ryHyr, template=pao%matrix_BFGS, matrix_type="N")
     306         276 :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ryH, matrix_yr, 0.0_dp, matrix_ryHyr)
     307             : 
     308             :          ! rr = MATMUL(r,TRANSPOSE(r))
     309         276 :          CALL dbcsr_create(matrix_rr, template=pao%matrix_BFGS, matrix_type="N")
     310         276 :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_r, matrix_r, 0.0_dp, matrix_rr)
     311             : 
     312             :          ! H = H - Hyr/ry - ryH/ry + ryHyr/(ry**2) + rr/ry
     313         276 :          CALL dbcsr_add(pao%matrix_BFGS, matrix_HYR, 1.0_dp, -1.0_dp/trace_ry)
     314         276 :          CALL dbcsr_add(pao%matrix_BFGS, matrix_ryH, 1.0_dp, -1.0_dp/trace_ry)
     315         276 :          CALL dbcsr_add(pao%matrix_BFGS, matrix_ryHyr, 1.0_dp, +1.0_dp/(trace_ry**2))
     316         276 :          CALL dbcsr_add(pao%matrix_BFGS, matrix_rr, 1.0_dp, +1.0_dp/trace_ry)
     317             : 
     318             :          ! clean up
     319         276 :          CALL dbcsr_release(matrix_y)
     320         276 :          CALL dbcsr_release(matrix_s)
     321         276 :          CALL dbcsr_release(matrix_r)
     322         276 :          CALL dbcsr_release(matrix_Hy)
     323         276 :          CALL dbcsr_release(matrix_yr)
     324         276 :          CALL dbcsr_release(matrix_Hyr)
     325         276 :          CALL dbcsr_release(matrix_ryH)
     326         276 :          CALL dbcsr_release(matrix_ryHyr)
     327         276 :          CALL dbcsr_release(matrix_rr)
     328             :       END IF
     329             : 
     330             :       ! approximate condition of Hessian
     331             :       !TODO: good setting for arnoldi?
     332             :       CALL arnoldi_extremal(pao%matrix_BFGS, eval_max, eval_min, max_iter=100, &
     333         288 :                             threshold=1e-2_dp, converged=arnoldi_converged)
     334         288 :       IF (arnoldi_converged) THEN
     335         432 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| evals of inv. Hessian: min, max, max/min", &
     336         288 :             eval_min, eval_max, eval_max/eval_min
     337             :       ELSE
     338           0 :          IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| arnoldi of inv. Hessian did not converged."
     339             :       END IF
     340             : 
     341             :       ! calculate new direction
     342             :       ! d = MATMUL(H, -g)
     343             :       CALL dbcsr_multiply("N", "N", -1.0_dp, pao%matrix_BFGS, matrix_G, &
     344         288 :                           0.0_dp, matrix_D, retain_sparsity=.TRUE.)
     345             : 
     346         288 :       CALL timestop(handle)
     347         288 :    END SUBROUTINE pao_opt_newdir_bfgs
     348             : 
     349             : END MODULE pao_optimizer