(git:ed6f26b)
Loading...
Searching...
No Matches
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-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! **************************************************************************************************
14MODULE 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
51CONTAINS
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_dp)
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
331END MODULE almo_scf_lbfgs_types
332
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
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
almo_scf_lbfgs_types::lbfgs_history_type
Definition almo_scf_lbfgs_types.F:41
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174