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