(git:ccc2433)
almo_scf_lbfgs_types.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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_log_handling, ONLY: cp_get_default_logger,&
17  cp_logger_get_default_unit_nr
18  !USE cp_log_handling, ONLY: cp_to_string
19  USE dbcsr_api, ONLY: dbcsr_add,&
20  dbcsr_copy,&
21  dbcsr_create,&
22  dbcsr_dot,&
23  dbcsr_release,&
24  dbcsr_scale,&
25  dbcsr_type
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  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  CALL lbfgs_history_push(history, variable, vartype=1, action=1)
65  CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
66 
67  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  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  CALL lbfgs_history_push(history, variable, vartype=1, action=2)
86  CALL lbfgs_history_push(history, gradient, vartype=2, action=2)
87  ! compute rho for the last stored value
88  CALL lbfgs_history_last_rho(history)
89 
90  CALL lbfgs_history_direction(history, gradient, direction)
91 
92  ! action 1 will seed the next storage cell
93  CALL lbfgs_history_push(history, variable, vartype=1, action=1)
94  CALL lbfgs_history_push(history, gradient, vartype=2, action=1)
95 
96  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  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  nallocate = max(1, nstore)
112  history%nstore = nallocate
113  history%istore(:) = 0 ! total number of action-2 pushes
114  ALLOCATE (history%matrix(nspins, nallocate, 2))
115  ALLOCATE (history%rho(nspins, nallocate))
116 
117  END SUBROUTINE lbfgs_create
118 
119 ! **************************************************************************************************
120 !> \brief release the bfgs history
121 !> \param history ...
122 ! **************************************************************************************************
123  SUBROUTINE lbfgs_release(history)
124  TYPE(lbfgs_history_type), INTENT(INOUT) :: history
125 
126  INTEGER :: ispin, istore, ivartype
127 
128  ! delete history
129  DO ispin = 1, SIZE(history%matrix, 1)
130  DO ivartype = 1, 2
131  DO istore = 1, min(history%istore(ivartype) + 1, history%nstore)
132  !WRITE(*,*) "ZREL: ispin,istore,vartype", ispin, istore, ivartype
133  CALL dbcsr_release(history%matrix(ispin, istore, ivartype))
134  END DO
135  END DO
136  END DO
137  DEALLOCATE (history%matrix)
138  DEALLOCATE (history%rho)
139 
140  END SUBROUTINE lbfgs_release
141 
142 ! **************************************************************************************************
143 !> \brief once all data in the last cell is stored, compute rho
144 !> \param history ...
145 ! **************************************************************************************************
146  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  DO ispin = 1, SIZE(history%matrix, 1)
160 
161  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  history%rho(ispin, istore))
165 
166  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  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  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  history%istore(vartype) = history%istore(vartype) + 1
203 
204  DO ispin = 1, SIZE(history%matrix, 1)
205 
206  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  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  template=matrix(ispin))
216  !IF (unit_nr > 0) THEN
217  ! WRITE (unit_nr, *) "Creating new matrix..."
218  !END IF
219  END IF
220 
221  IF (action .EQ. 1) THEN
222  CALL dbcsr_copy(history%matrix(ispin, istore, vartype), matrix(ispin))
223  ELSE
224  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  IF (action .EQ. 1) THEN
232  history%istore(vartype) = history%istore(vartype) - 1
233  END IF
234 
235  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  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  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  IF (history%istore(1) .NE. history%istore(2)) THEN
262  cpabort("BFGS APIs are not used correctly")
263  END IF
264 
265  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  ALLOCATE (alpha(nterms))
271 
272  DO ispin = 1, SIZE(history%matrix, 1)
273 
274  CALL dbcsr_create(q, template=gradient(ispin))
275 
276  CALL dbcsr_copy(q, gradient(ispin))
277 
278  ! loop over all stored items
279  DO iterm = 1, nterms
280 
281  ! location: from recent to oldest stored
282  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  CALL dbcsr_dot(history%matrix(ispin, istore, 1), q, alpha(iterm))
289  alpha(iterm) = history%rho(ispin, istore)*alpha(iterm)
290  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  IF (iterm .EQ. 1) THEN
295  CALL dbcsr_dot(history%matrix(ispin, istore, 2), history%matrix(ispin, istore, 2), gammak)
296  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  CALL dbcsr_scale(q, gammak)
306 
307  ! loop over all stored items
308  DO iterm = nterms, 1, -1
309 
310  ! location: from oldest to recent stored
311  istore = mod(history%istore(1) - iterm, history%nstore) + 1
312 
313  CALL dbcsr_dot(history%matrix(ispin, istore, 2), q, beta)
314  beta = history%rho(ispin, istore)*beta
315  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  CALL dbcsr_scale(q, -1.0)
321  CALL dbcsr_copy(direction(ispin), q)
322 
323  CALL dbcsr_release(q)
324 
325  END DO !ispin
326 
327  DEALLOCATE (alpha)
328 
329  END SUBROUTINE lbfgs_history_direction
330 
331 END MODULE almo_scf_lbfgs_types
332 
almo_scf_lbfgs_types
Limited memory BFGS.
Definition: almo_scf_lbfgs_types.F:14
almo_scf_lbfgs_types::lbfgs_create
subroutine, public lbfgs_create(history, nspins, nstore)
create history storage for limited memory bfgs
Definition: almo_scf_lbfgs_types.F:105
almo_scf_lbfgs_types::lbfgs_seed
subroutine, public lbfgs_seed(history, variable, gradient)
interface subroutine to store the first variable/gradient pair
Definition: almo_scf_lbfgs_types.F:60
almo_scf_lbfgs_types::lbfgs_release
subroutine, public lbfgs_release(history)
release the bfgs history
Definition: almo_scf_lbfgs_types.F:124
almo_scf_lbfgs_types::lbfgs_get_direction
subroutine, public lbfgs_get_direction(history, variable, gradient, direction)
interface subroutine to store a variable/gradient pair and predict direction
Definition: almo_scf_lbfgs_types.F:79
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition: cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34