Line data Source code
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 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)
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 :
|
