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