LCOV - code coverage report
Current view: top level - src - almo_scf_lbfgs_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.3 % 74 72
Test Date: 2025-07-25 12:55:17 Functions: 77.8 % 9 7

            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 Limited memory BFGS
      10              : !> \par History
      11              : !>       2019.10 created [Rustam Z Khaliullin]
      12              : !> \author Rustam Z Khaliullin
      13              : ! **************************************************************************************************
      14              : MODULE almo_scf_lbfgs_types
      15              :    !USE cp_external_control,             ONLY: external_control
      16              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      17              :                                               dbcsr_copy,&
      18              :                                               dbcsr_create,&
      19              :                                               dbcsr_release,&
      20              :                                               dbcsr_scale,&
      21              :                                               dbcsr_type
      22              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot
      23              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24              :                                               cp_logger_get_default_unit_nr
      25              :    !USE cp_log_handling,                 ONLY: cp_to_string
      26              :    USE kinds,                           ONLY: dp
      27              : #include "./base/base_uses.f90"
      28              : 
      29              :    IMPLICIT NONE
      30              : 
      31              :    PRIVATE
      32              : 
      33              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_lbfgs_types'
      34              : 
      35              :    PUBLIC :: lbfgs_seed, &
      36              :              lbfgs_create, &
      37              :              lbfgs_release, &
      38              :              lbfgs_get_direction, &
      39              :              lbfgs_history_type
      40              : 
      41              :    TYPE lbfgs_history_type
      42              :       INTEGER :: nstore = 0
      43              :       ! istore counts the total number of action=2 pushes
      44              :       ! istore is designed to become more than nstore eventually
      45              :       ! there are two counters: the main variable and gradient
      46              :       INTEGER, DIMENSION(2) :: istore = 0
      47              :       TYPE(dbcsr_type), DIMENSION(:, :, :), ALLOCATABLE :: matrix
      48              :       REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rho
      49              :    END TYPE lbfgs_history_type
      50              : 
      51              : CONTAINS
      52              : 
      53              : ! **************************************************************************************************
      54              : !> \brief interface subroutine to store the first variable/gradient pair
      55              : !> \param history ...
      56              : !> \param variable ...
      57              : !> \param gradient ...
      58              : ! **************************************************************************************************
      59            6 :    SUBROUTINE lbfgs_seed(history, variable, gradient)
      60              : 
      61              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
      62              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: variable, gradient
      63              : 
      64            6 :       CALL lbfgs_history_push(history, variable, vartype=1, action=1)
      65            6 :       CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
      66              : 
      67            6 :    END SUBROUTINE lbfgs_seed
      68              : 
      69              : ! **************************************************************************************************
      70              : !> \brief interface subroutine to store a variable/gradient pair
      71              : !>        and predict direction
      72              : !> \param history ...
      73              : !> \param variable ...
      74              : !> \param gradient ...
      75              : !> \param direction ...
      76              : ! **************************************************************************************************
      77              : 
      78           24 :    SUBROUTINE lbfgs_get_direction(history, variable, gradient, direction)
      79              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
      80              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: variable, gradient
      81              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT)      :: direction
      82              : 
      83              :       ! action 2 will calculate delta = (new - old)
      84              :       ! in the last used storage cell
      85           24 :       CALL lbfgs_history_push(history, variable, vartype=1, action=2)
      86           24 :       CALL lbfgs_history_push(history, gradient, vartype=2, action=2)
      87              :       ! compute rho for the last stored value
      88           24 :       CALL lbfgs_history_last_rho(history)
      89              : 
      90           24 :       CALL lbfgs_history_direction(history, gradient, direction)
      91              : 
      92              :       ! action 1 will seed the next storage cell
      93           24 :       CALL lbfgs_history_push(history, variable, vartype=1, action=1)
      94           24 :       CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
      95              : 
      96           24 :    END SUBROUTINE lbfgs_get_direction
      97              : 
      98              : ! **************************************************************************************************
      99              : !> \brief create history storage for limited memory bfgs
     100              : !> \param history ...
     101              : !> \param nspins ...
     102              : !> \param nstore ...
     103              : ! **************************************************************************************************
     104            8 :    SUBROUTINE lbfgs_create(history, nspins, nstore)
     105              : 
     106              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
     107              :       INTEGER, INTENT(IN)                                :: nspins, nstore
     108              : 
     109              :       INTEGER                                            :: nallocate
     110              : 
     111            8 :       nallocate = MAX(1, nstore)
     112            8 :       history%nstore = nallocate
     113           24 :       history%istore(:) = 0  ! total number of action-2 pushes
     114          376 :       ALLOCATE (history%matrix(nspins, nallocate, 2))
     115           32 :       ALLOCATE (history%rho(nspins, nallocate))
     116              : 
     117            8 :    END SUBROUTINE lbfgs_create
     118              : 
     119              : ! **************************************************************************************************
     120              : !> \brief release the bfgs history
     121              : !> \param history ...
     122              : ! **************************************************************************************************
     123            8 :    SUBROUTINE lbfgs_release(history)
     124              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
     125              : 
     126              :       INTEGER                                            :: ispin, istore, ivartype
     127              : 
     128              :       ! delete history
     129           16 :       DO ispin = 1, SIZE(history%matrix, 1)
     130           32 :          DO ivartype = 1, 2
     131           88 :             DO istore = 1, MIN(history%istore(ivartype) + 1, history%nstore)
     132              :                !WRITE(*,*) "ZREL: ispin,istore,vartype", ispin, istore, ivartype
     133           80 :                CALL dbcsr_release(history%matrix(ispin, istore, ivartype))
     134              :             END DO
     135              :          END DO
     136              :       END DO
     137            8 :       DEALLOCATE (history%matrix)
     138            8 :       DEALLOCATE (history%rho)
     139              : 
     140            8 :    END SUBROUTINE lbfgs_release
     141              : 
     142              : ! **************************************************************************************************
     143              : !> \brief once all data in the last cell is stored, compute rho
     144              : !> \param history ...
     145              : ! **************************************************************************************************
     146           24 :    SUBROUTINE lbfgs_history_last_rho(history)
     147              : 
     148              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
     149              : 
     150              :       INTEGER                                            :: ispin, istore
     151              : 
     152              :       !logger => cp_get_default_logger()
     153              :       !IF (logger%para_env%is_source()) THEN
     154              :       !   unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     155              :       !ELSE
     156              :       !   unit_nr = -1
     157              :       !ENDIF
     158              : 
     159           48 :       DO ispin = 1, SIZE(history%matrix, 1)
     160              : 
     161           24 :          istore = MOD(history%istore(1) - 1, history%nstore) + 1
     162              :          CALL dbcsr_dot(history%matrix(ispin, istore, 1), &
     163              :                         history%matrix(ispin, istore, 2), &
     164           24 :                         history%rho(ispin, istore))
     165              : 
     166           48 :          history%rho(ispin, istore) = 1.0_dp/history%rho(ispin, istore)
     167              : 
     168              :          !IF (unit_nr > 0) THEN
     169              :          !   WRITE (unit_nr, *) "Rho in cell ", istore, " is computed ", history%rho(ispin, istore)
     170              :          !ENDIF
     171              : 
     172              :       END DO ! ispin
     173              : 
     174           24 :    END SUBROUTINE lbfgs_history_last_rho
     175              : 
     176              : ! **************************************************************************************************
     177              : !> \brief store data in history
     178              : !>  vartype - which data piece to store: 1 - variable, 2 - gradient
     179              : !>  operation - what to do: 1 - erase existing and store new
     180              : !>                          2 - store = new - existing
     181              : !> \param history ...
     182              : !> \param matrix ...
     183              : !> \param vartype ...
     184              : !> \param action ...
     185              : ! **************************************************************************************************
     186          108 :    SUBROUTINE lbfgs_history_push(history, matrix, vartype, action)
     187              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
     188              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: matrix
     189              :       INTEGER, INTENT(IN)                                :: vartype, action
     190              : 
     191              :       INTEGER                                            :: ispin, istore
     192              : 
     193              :       !logger => cp_get_default_logger()
     194              :       !IF (logger%para_env%is_source()) THEN
     195              :       !   unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     196              :       !ELSE
     197              :       !   unit_nr = -1
     198              :       !ENDIF
     199              : 
     200              :       ! increase the counter: it moves the pointer to the next cell
     201              :       ! for action==1 this is a "pretend" increase; the pointer will be moved back in the end
     202          108 :       history%istore(vartype) = history%istore(vartype) + 1
     203              : 
     204          216 :       DO ispin = 1, SIZE(history%matrix, 1)
     205              : 
     206          108 :          istore = MOD(history%istore(vartype) - 1, history%nstore) + 1
     207              :          !IF (unit_nr > 0) THEN
     208              :          !   WRITE (unit_nr, *) "Action ", action, " modifying cell ", istore
     209              :          !END IF
     210              : 
     211          108 :          IF (history%istore(vartype) <= history%nstore .AND. &
     212              :              action .EQ. 1) THEN
     213              :             !WRITE(*,*) "ZCRE: ispin,istore,vartype", ispin, istore, vartype
     214              :             CALL dbcsr_create(history%matrix(ispin, istore, vartype), &
     215           60 :                               template=matrix(ispin))
     216              :             !IF (unit_nr > 0) THEN
     217              :             !   WRITE (unit_nr, *) "Creating new matrix..."
     218              :             !END IF
     219              :          END IF
     220              : 
     221          216 :          IF (action .EQ. 1) THEN
     222           60 :             CALL dbcsr_copy(history%matrix(ispin, istore, vartype), matrix(ispin))
     223              :          ELSE
     224           48 :             CALL dbcsr_add(history%matrix(ispin, istore, vartype), matrix(ispin), -1.0_dp, 1.0_dp)
     225              :          END IF
     226              : 
     227              :       END DO ! ispin
     228              : 
     229              :       ! allow the pointer to move forward only if deltas are stored (action==2)
     230              :       ! otherwise return the pointer to the previous value
     231          108 :       IF (action .EQ. 1) THEN
     232           60 :          history%istore(vartype) = history%istore(vartype) - 1
     233              :       END IF
     234              : 
     235          108 :    END SUBROUTINE lbfgs_history_push
     236              : 
     237              : ! **************************************************************************************************
     238              : !> \brief use history data to construct dir = -Hinv.grad
     239              : !> \param history ...
     240              : !> \param gradient ...
     241              : !> \param direction ...
     242              : ! **************************************************************************************************
     243           24 :    SUBROUTINE lbfgs_history_direction(history, gradient, direction)
     244              : 
     245              :       TYPE(lbfgs_history_type), INTENT(INOUT)            :: history
     246              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: gradient
     247              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT)      :: direction
     248              : 
     249              :       INTEGER                                            :: ispin, istore, iterm, nterms
     250              :       REAL(KIND=dp)                                      :: beta, gammak
     251           24 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: alpha
     252              :       TYPE(dbcsr_type)                                   :: q
     253              : 
     254              :       !logger => cp_get_default_logger()
     255              :       !IF (logger%para_env%is_source()) THEN
     256              :       !   unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     257              :       !ELSE
     258              :       !   unit_nr = -1
     259              :       !ENDIF
     260              : 
     261           24 :       IF (history%istore(1) .NE. history%istore(2)) THEN
     262            0 :          CPABORT("BFGS APIs are not used correctly")
     263              :       END IF
     264              : 
     265           24 :       nterms = MIN(history%istore(1), history%nstore)
     266              :       !IF (unit_nr > 0) THEN
     267              :       !   WRITE (unit_nr, *) "L-BFGS terms used: ", nterms
     268              :       !END IF
     269              : 
     270           72 :       ALLOCATE (alpha(nterms))
     271              : 
     272           48 :       DO ispin = 1, SIZE(history%matrix, 1)
     273              : 
     274           24 :          CALL dbcsr_create(q, template=gradient(ispin))
     275              : 
     276           24 :          CALL dbcsr_copy(q, gradient(ispin))
     277              : 
     278              :          ! loop over all stored items
     279          108 :          DO iterm = 1, nterms
     280              : 
     281              :             ! location: from recent to oldest stored
     282           84 :             istore = MOD(history%istore(1) - iterm, history%nstore) + 1
     283              : 
     284              :             !IF (unit_nr > 0) THEN
     285              :             !   WRITE (unit_nr, *) "Record locator: ", istore
     286              :             !END IF
     287              : 
     288           84 :             CALL dbcsr_dot(history%matrix(ispin, istore, 1), q, alpha(iterm))
     289           84 :             alpha(iterm) = history%rho(ispin, istore)*alpha(iterm)
     290           84 :             CALL dbcsr_add(q, history%matrix(ispin, istore, 2), 1.0_dp, -alpha(iterm))
     291              : 
     292              :             ! use the most recent term to
     293              :             ! compute gamma_k, Nocedal (7.20) and then get H0
     294          108 :             IF (iterm .EQ. 1) THEN
     295           24 :                CALL dbcsr_dot(history%matrix(ispin, istore, 2), history%matrix(ispin, istore, 2), gammak)
     296           24 :                gammak = 1.0_dp/(gammak*history%rho(ispin, istore))
     297              :                !IF (unit_nr > 0) THEN
     298              :                !   WRITE (unit_nr, *) "Gamma_k: ", gammak
     299              :                !END IF
     300              :             END IF
     301              : 
     302              :          END DO ! iterm, first loop from recent to oldest
     303              : 
     304              :          ! now q stores Nocedal's r = (gamma_k*I).q
     305           24 :          CALL dbcsr_scale(q, gammak)
     306              : 
     307              :          ! loop over all stored items
     308          108 :          DO iterm = nterms, 1, -1
     309              : 
     310              :             ! location: from oldest to recent stored
     311           84 :             istore = MOD(history%istore(1) - iterm, history%nstore) + 1
     312              : 
     313           84 :             CALL dbcsr_dot(history%matrix(ispin, istore, 2), q, beta)
     314           84 :             beta = history%rho(ispin, istore)*beta
     315          108 :             CALL dbcsr_add(q, history%matrix(ispin, istore, 1), 1.0_dp, alpha(iterm) - beta)
     316              : 
     317              :          END DO ! iterm, forst loop from recent to oldest
     318              : 
     319              :          !RZK-warning: unclear whether q should be multiplied by minus one
     320           24 :          CALL dbcsr_scale(q, -1.0_dp)
     321           24 :          CALL dbcsr_copy(direction(ispin), q)
     322              : 
     323           48 :          CALL dbcsr_release(q)
     324              : 
     325              :       END DO !ispin
     326              : 
     327           24 :       DEALLOCATE (alpha)
     328              : 
     329           24 :    END SUBROUTINE lbfgs_history_direction
     330              : 
     331            0 : END MODULE almo_scf_lbfgs_types
     332              : 
        

Generated by: LCOV version 2.0-1