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