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 Routines for Geometry optimization using BFGS algorithm
10 : !> \par History
11 : !> Module modified by Pierre-André Cazade [pcazade] 01.2020 - University of Limerick.
12 : !> Modifications enable Space Group Symmetry.
13 : ! **************************************************************************************************
14 : MODULE bfgs_optimizer
15 :
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 : USE atomic_kind_types, ONLY: get_atomic_kind,&
18 : get_atomic_kind_set
19 : USE bibliography, ONLY: Lindh1995,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type,&
22 : pbc
23 : USE constraint_fxd, ONLY: fix_atom_control
24 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
25 : cp_blacs_env_release,&
26 : cp_blacs_env_type
27 : USE cp_external_control, ONLY: external_control
28 : USE cp_files, ONLY: close_file,&
29 : open_file
30 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
31 : cp_fm_matvec,&
32 : cp_fm_transpose
33 : USE cp_fm_diag, ONLY: choose_eigv_solver
34 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 : cp_fm_struct_release,&
36 : cp_fm_struct_type
37 : USE cp_fm_types, ONLY: cp_fm_create,&
38 : cp_fm_get_info,&
39 : cp_fm_read_unformatted,&
40 : cp_fm_release,&
41 : cp_fm_set_all,&
42 : cp_fm_to_fm,&
43 : cp_fm_type,&
44 : cp_fm_write_unformatted
45 : USE cp_log_handling, ONLY: cp_get_default_logger,&
46 : cp_logger_type,&
47 : cp_to_string
48 : USE cp_output_handling, ONLY: cp_iterate,&
49 : cp_p_file,&
50 : cp_print_key_finished_output,&
51 : cp_print_key_should_output,&
52 : cp_print_key_unit_nr
53 : USE cp_subsys_types, ONLY: cp_subsys_get,&
54 : cp_subsys_type
55 : USE force_env_types, ONLY: force_env_get,&
56 : force_env_type
57 : USE global_types, ONLY: global_environment_type
58 : USE gopt_f_methods, ONLY: &
59 : cp_eval_at, gopt_f_ii, gopt_f_io, gopt_f_io_finalize, gopt_f_io_init, &
60 : gopt_f_project_coord, gopt_f_project_step, gopt_f_project_vector, gopt_f_space_group_init, &
61 : gopt_f_space_group_show, print_geo_opt_header, print_geo_opt_nc
62 : USE gopt_f_types, ONLY: gopt_f_type
63 : USE gopt_param_types, ONLY: gopt_param_type
64 : USE ieee_arithmetic, ONLY: ieee_is_finite
65 : USE input_constants, ONLY: bfgs_direct,&
66 : bfgs_inverse,&
67 : default_cell_method_id,&
68 : default_ts_method_id
69 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
70 : section_vals_type,&
71 : section_vals_val_get,&
72 : section_vals_val_set
73 : USE iso_fortran_env, ONLY: iostat_end
74 : USE kinds, ONLY: default_path_length,&
75 : dp
76 : USE machine, ONLY: m_flush,&
77 : m_walltime
78 : USE mathlib, ONLY: det_3x3,&
79 : inv_3x3
80 : USE message_passing, ONLY: mp_para_env_type
81 : USE parallel_gemm_api, ONLY: parallel_gemm
82 : USE particle_list_types, ONLY: particle_list_type
83 : #include "../base/base_uses.f90"
84 :
85 : IMPLICIT NONE
86 : PRIVATE
87 :
88 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bfgs_optimizer'
89 : ! Parameters of the Lindh 1995 model Hessian (atom rows H/C/O-N, pair
90 : ! columns). Screening constants for construct_initial_hess.
91 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_alpha(3, 3) = &
92 : RESHAPE([1.0_dp, 0.3494_dp, 0.3494_dp, &
93 : 0.3949_dp, 0.2800_dp, 0.2800_dp, &
94 : 0.3949_dp, 0.2800_dp, 0.1800_dp], [3, 3])
95 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_r0(3, 3) = &
96 : RESHAPE([1.35_dp, 2.10_dp, 2.53_dp, &
97 : 2.10_dp, 2.87_dp, 3.40_dp, &
98 : 2.53_dp, 3.40_dp, 3.40_dp], [3, 3])
99 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_cutoff = 6.0_dp
100 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_stretch_factor = 0.45_dp
101 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_angle_factor = 0.15_dp
102 : REAL(KIND=dp), PARAMETER, PRIVATE :: lindh_denom_floor = 0.01_dp
103 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
104 :
105 : ! Both matrix representations share one restart format and the existing file name.
106 : CHARACTER(LEN=16), PARAMETER :: bfgs_restart_magic = "CP2K BFGS STATE "
107 : INTEGER, PARAMETER :: bfgs_restart_version = 1
108 : REAL(KIND=dp), PARAMETER :: min_hessian_eigenvalue = 1.0e-4_dp
109 :
110 : PUBLIC :: bfgs_read_restart_header, geoopt_bfgs
111 :
112 : CONTAINS
113 :
114 : ! **************************************************************************************************
115 : !> \brief Main driver for BFGS geometry optimizations
116 : !> \param force_env ...
117 : !> \param gopt_param ...
118 : !> \param globenv ...
119 : !> \param geo_section ...
120 : !> \param gopt_env ...
121 : !> \param x0 ...
122 : !> \par History
123 : !> 01.2020 modified to perform Space Group Symmetry [pcazade]
124 : ! **************************************************************************************************
125 905 : RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
126 :
127 : TYPE(force_env_type), POINTER :: force_env
128 : TYPE(gopt_param_type), POINTER :: gopt_param
129 : TYPE(global_environment_type), POINTER :: globenv
130 : TYPE(section_vals_type), POINTER :: geo_section
131 : TYPE(gopt_f_type), POINTER :: gopt_env
132 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
133 :
134 : CHARACTER(len=*), PARAMETER :: routineN = 'geoopt_bfgs'
135 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
136 :
137 : CHARACTER(LEN=5) :: wildcard
138 : INTEGER :: handle, info, iter_nr, its, &
139 : matrix_method, maxiter, ndf, &
140 : output_unit
141 : LOGICAL :: conv, hesrest, hess_reset, inverse, &
142 : shell_present, should_stop, &
143 : use_mod_hes, use_rfo
144 : REAL(KIND=dp) :: ediff, emin, eold, etot, pred, rad, rat, &
145 : step, t_diff, t_now, t_old
146 905 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
147 905 : REAL(KIND=dp), DIMENSION(:), POINTER :: g
148 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
149 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
150 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_hes
151 : TYPE(cp_fm_type) :: eigvec_mat, hess_mat, hess_tmp
152 : TYPE(cp_logger_type), POINTER :: logger
153 : TYPE(cp_subsys_type), POINTER :: subsys
154 : TYPE(mp_para_env_type), POINTER :: para_env
155 : TYPE(section_vals_type), POINTER :: print_key, root_section
156 :
157 905 : NULLIFY (logger, g, blacs_env)
158 1810 : logger => cp_get_default_logger()
159 905 : para_env => force_env%para_env
160 905 : root_section => force_env%root_section
161 905 : t_old = m_walltime()
162 :
163 905 : CALL timeset(routineN, handle)
164 905 : CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
165 905 : print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
166 905 : CALL section_vals_val_get(geo_section, "BFGS%METHOD", i_val=matrix_method)
167 905 : inverse = matrix_method == bfgs_inverse
168 905 : maxiter = gopt_param%max_iter
169 905 : conv = .FALSE.
170 905 : rat = 0.0_dp
171 905 : wildcard = " BFGS"
172 :
173 : ! Stop if not yet implemented
174 905 : SELECT CASE (gopt_env%type_id)
175 : CASE (default_ts_method_id)
176 905 : CPABORT("BFGS method not yet working with DIMER")
177 : END SELECT
178 :
179 905 : CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
180 905 : CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
181 905 : CALL section_vals_val_get(geo_section, "BFGS%RESTART_HESSIAN", l_val=hesrest)
182 905 : IF (inverse .AND. use_rfo) THEN
183 : CALL cp_abort(__LOCATION__, &
184 0 : "BFGS METHOD INVERSE does not support USE_RAT_FUN_OPT")
185 : END IF
186 : IF (inverse) THEN
187 10 : IF (.NOT. ieee_is_finite(rad) .OR. rad <= 0.0_dp) THEN
188 : CALL cp_abort(__LOCATION__, &
189 0 : "BFGS TRUST_RADIUS must be finite and positive")
190 : END IF
191 : END IF
192 : output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
193 905 : extension=".geoLog")
194 905 : IF (output_unit > 0) THEN
195 476 : IF (inverse) THEN
196 5 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "BFGS| Matrix update method: INVERSE"
197 : ELSE
198 471 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "BFGS| Matrix update method: DIRECT"
199 : END IF
200 476 : IF (use_rfo) THEN
201 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
202 6 : "BFGS| Use rational function optimization for step estimation: ", "YES"
203 : ELSE
204 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
205 470 : "BFGS| Use rational function optimization for step estimation: ", " NO"
206 : END IF
207 476 : IF (use_mod_hes) THEN
208 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
209 410 : "BFGS| Use model Hessian for initial guess: ", "YES"
210 : ELSE
211 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
212 66 : "BFGS| Use model Hessian for initial guess: ", " NO"
213 : END IF
214 476 : IF (hesrest) THEN
215 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
216 4 : "BFGS| Restart Hessian: ", "YES"
217 : ELSE
218 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
219 472 : "BFGS| Restart Hessian: ", " NO"
220 : END IF
221 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
222 476 : "BFGS| Trust radius: ", rad
223 : END IF
224 :
225 905 : ndf = SIZE(x0)
226 905 : IF (ndf > 3000 .AND. .NOT. inverse) THEN
227 : CALL cp_warn(__LOCATION__, &
228 : "The dimension of the Hessian matrix ("// &
229 : TRIM(ADJUSTL(cp_to_string(ndf)))//") is greater than 3000. "// &
230 : "The diagonalisation of the full Hessian matrix needed for BFGS "// &
231 : "is computationally expensive. You should consider to use the linear "// &
232 0 : "scaling variant L-BFGS instead.")
233 : END IF
234 :
235 : ! Initialize hessian (hes = unitary matrix or model hessian )
236 : CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
237 905 : globenv%blacs_repeatable)
238 : CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
239 905 : nrow_global=ndf, ncol_global=ndf)
240 905 : CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
241 :
242 905 : CALL force_env_get(force_env=force_env, subsys=subsys)
243 905 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
244 905 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
245 905 : IF (use_mod_hes) THEN
246 773 : IF (shell_present) THEN
247 : CALL cp_warn(__LOCATION__, &
248 : "No model Hessian is available for core-shell models. "// &
249 4 : "A unit matrix is used as the initial Hessian.")
250 4 : use_mod_hes = .FALSE.
251 : END IF
252 773 : IF (gopt_env%type_id == default_cell_method_id) THEN
253 : CALL cp_warn(__LOCATION__, &
254 : "No model Hessian is available for cell optimizations. "// &
255 0 : "A unit matrix is used as the initial Hessian.")
256 0 : use_mod_hes = .FALSE.
257 : END IF
258 : END IF
259 :
260 : ! hess_mat stores B for DIRECT and H = B^{-1} for INVERSE.
261 : ! DIRECT keeps the spectral work arrays for step construction; INVERSE
262 : ! only needs them while building a model-Hessian initial guess.
263 905 : IF (.NOT. inverse .OR. (use_mod_hes .AND. .NOT. hesrest)) THEN
264 901 : CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
265 901 : CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
266 2703 : ALLOCATE (eigval(ndf))
267 1802 : eigval(:) = zero
268 : END IF
269 :
270 905 : IF (inverse) THEN
271 10 : IF (use_mod_hes .AND. .NOT. hesrest) THEN
272 : CALL build_initial_hess(hess_mat, hess_tmp, eigvec_mat, eigval, force_env, &
273 6 : use_mod_hes, .TRUE.)
274 : ELSE
275 4 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
276 : END IF
277 : ELSE
278 : CALL build_initial_hess(hess_mat, hess_tmp, eigvec_mat, eigval, force_env, &
279 895 : use_mod_hes, .FALSE.)
280 : END IF
281 :
282 905 : IF (inverse .AND. ALLOCATED(eigval)) THEN
283 6 : CALL cp_fm_release(eigvec_mat)
284 6 : CALL cp_fm_release(hess_tmp)
285 6 : DEALLOCATE (eigval)
286 : END IF
287 :
288 : ! Restore the reference coordinates before symmetry setup and the first evaluation.
289 905 : IF (hesrest) CALL read_bfgs_hessian(geo_section, hess_mat, logger, gopt_env, x0, matrix_method)
290 :
291 2715 : ALLOCATE (xold(ndf))
292 25493 : xold(:) = x0(:)
293 :
294 1810 : ALLOCATE (g(ndf))
295 25493 : g(:) = 0.0_dp
296 :
297 1810 : ALLOCATE (gold(ndf))
298 905 : gold(:) = 0.0_dp
299 :
300 1810 : ALLOCATE (dx(ndf))
301 905 : dx(:) = 0.0_dp
302 :
303 1810 : ALLOCATE (dg(ndf))
304 905 : dg(:) = 0.0_dp
305 :
306 1810 : ALLOCATE (work(ndf))
307 905 : work(:) = 0.0_dp
308 :
309 1810 : ALLOCATE (dr(ndf))
310 905 : dr(:) = 0.0_dp
311 :
312 : ! find space_group
313 905 : CALL gopt_f_space_group_init(gopt_env, x0, geo_section, output_unit)
314 :
315 : ! Geometry optimization starts now
316 905 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
317 905 : CALL print_geo_opt_header(gopt_env, output_unit, wildcard)
318 :
319 : ! Calculate Energy & Gradients
320 : CALL cp_eval_at(gopt_env, x0, etot, g, master=gopt_env%force_env%para_env%mepos, &
321 905 : para_env=gopt_env%force_env%para_env)
322 :
323 : ! Symmetrize coordinates and forces
324 905 : CALL gopt_f_project_coord(gopt_env, x0)
325 905 : CALL gopt_f_project_vector(gopt_env, g)
326 :
327 : ! Print info at time 0
328 905 : emin = etot
329 905 : t_now = m_walltime()
330 905 : t_diff = t_now - t_old
331 905 : t_old = t_now
332 905 : CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
333 4392 : DO its = iter_nr + 1, maxiter
334 4382 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
335 4382 : CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
336 4382 : CALL gopt_f_ii(its, output_unit)
337 :
338 4382 : IF ((its - iter_nr) > 1) THEN
339 : ! Symmetrize old coordinates and old forces.
340 3487 : CALL gopt_f_project_coord(gopt_env, xold)
341 3487 : CALL gopt_f_project_vector(gopt_env, gold)
342 :
343 : ! cp_eval_at may project the proposed geometry: use the actual displacement.
344 326215 : dx(:) = x0(:) - xold(:)
345 326215 : dg(:) = g(:) - gold(:)
346 3487 : IF (inverse) THEN
347 502 : CALL bfgs_inverse_update(ndf, dx, dg, hess_mat, work)
348 : ELSE
349 2985 : CALL bfgs(ndf, dx, dg, hess_mat, work, para_env, output_unit, hess_reset)
350 2985 : IF (hess_reset) THEN
351 : CALL build_initial_hess(hess_mat, hess_tmp, eigvec_mat, eigval, force_env, &
352 0 : use_mod_hes, .FALSE.)
353 : END IF
354 : END IF
355 :
356 3487 : CALL gopt_f_project_vector(gopt_env, dx)
357 3487 : CALL gopt_f_project_vector(gopt_env, dg)
358 :
359 3487 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
360 3162 : CALL write_bfgs_hessian(geo_section, hess_mat, logger, gopt_env, matrix_method)
361 : END IF
362 : END IF
363 :
364 : ! Symmetrize coordinates and forces
365 4382 : CALL gopt_f_project_coord(gopt_env, x0)
366 4382 : CALL gopt_f_project_vector(gopt_env, g)
367 :
368 : ! Setting the present positions & gradients as old
369 351554 : xold(:) = x0
370 351554 : gold(:) = g
371 :
372 4382 : IF (inverse) THEN
373 512 : CALL bfgs_inverse_step(hess_mat, g, dr)
374 : ELSE
375 : ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
376 3870 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
377 :
378 3870 : CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info)
379 :
380 : ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?)
381 3870 : IF (info /= 0) THEN
382 0 : IF (output_unit > 0) THEN
383 : WRITE (output_unit, *) &
384 0 : "BFGS: Matrix diagonalization failed, resetting Hessian to unity."
385 : END IF
386 0 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
387 0 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
388 0 : CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval)
389 : END IF
390 :
391 3870 : IF (use_rfo) THEN
392 530 : CALL set_hes_eig(ndf, eigval, work)
393 70634 : dx(:) = eigval
394 530 : CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
395 : END IF
396 3870 : CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
397 :
398 : END IF
399 4382 : CALL gopt_f_project_step(gopt_env, dr)
400 :
401 4382 : IF (inverse) THEN
402 : ! Avoid a zero-step division and initialize the reported scale on every step.
403 512 : rat = 1.0_dp
404 79712 : step = MAXVAL(ABS(dr))
405 512 : IF (step > rad) THEN
406 16 : rat = rad/step
407 1864 : dr(:) = rat*dr(:)
408 : step = rad
409 : END IF
410 : ELSE
411 3870 : CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
412 : END IF
413 :
414 : ! Update the atomic positions
415 351554 : x0 = x0 + dr
416 :
417 : ! Symmetrize coordinates
418 4382 : CALL gopt_f_project_coord(gopt_env, x0)
419 :
420 4382 : eold = etot
421 :
422 : ! Energy & Gradients at new step
423 : CALL cp_eval_at(gopt_env, x0, etot, g, master=gopt_env%force_env%para_env%mepos, &
424 4382 : para_env=gopt_env%force_env%para_env)
425 :
426 : ! Use the displacement actually accepted by constraints and cell projection.
427 351554 : dr(:) = x0(:) - xold(:)
428 351554 : step = MAXVAL(ABS(dr))
429 :
430 : ! This diagnostic requires B, not H. It is optional in gopt_f_io.
431 4382 : IF (.NOT. inverse) CALL energy_predict(ndf, work, hess_mat, dr, gold, pred, para_env)
432 :
433 4382 : ediff = etot - eold
434 :
435 : ! Symmetrize forces
436 4382 : CALL gopt_f_project_vector(gopt_env, g)
437 :
438 : ! check for an external exit command
439 4382 : CALL external_control(should_stop, "GEO", globenv=globenv)
440 4382 : IF (should_stop) EXIT
441 :
442 : ! Some IO and Convergence check
443 4382 : t_now = m_walltime()
444 4382 : t_diff = t_now - t_old
445 4382 : t_old = t_now
446 4382 : IF (inverse) THEN
447 : CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, &
448 : eold, emin, wildcard, gopt_param, ndf, dr, g, conv, &
449 512 : rat=rat, step=step, rad=rad, used_time=t_diff)
450 : ELSE
451 : CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, &
452 : eold, emin, wildcard, gopt_param, ndf, dr, g, conv, pred, rat, &
453 3870 : step, rad, used_time=t_diff)
454 : END IF
455 :
456 4382 : IF (conv .OR. (its == maxiter)) EXIT
457 3487 : IF (etot < emin) emin = etot
458 13156 : IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
459 : END DO
460 :
461 905 : IF (its == maxiter .AND. (.NOT. conv)) THEN
462 627 : CALL print_geo_opt_nc(gopt_env, output_unit)
463 : END IF
464 :
465 : ! show space_group
466 905 : CALL gopt_f_space_group_show(gopt_env, geo_section, output_unit)
467 :
468 : ! Write final information, if converged
469 905 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
470 905 : CALL write_bfgs_hessian(geo_section, hess_mat, logger, gopt_env, matrix_method)
471 : CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
472 905 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
473 :
474 905 : CALL cp_fm_struct_release(fm_struct_hes)
475 905 : CALL cp_fm_release(hess_mat)
476 905 : IF (.NOT. inverse) THEN
477 895 : CALL cp_fm_release(eigvec_mat)
478 895 : CALL cp_fm_release(hess_tmp)
479 : END IF
480 :
481 905 : CALL cp_blacs_env_release(blacs_env)
482 905 : DEALLOCATE (xold)
483 905 : DEALLOCATE (g)
484 905 : DEALLOCATE (gold)
485 905 : DEALLOCATE (dx)
486 905 : DEALLOCATE (dg)
487 905 : IF (ALLOCATED(eigval)) DEALLOCATE (eigval)
488 905 : DEALLOCATE (work)
489 905 : DEALLOCATE (dr)
490 :
491 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
492 905 : "PRINT%PROGRAM_RUN_INFO")
493 905 : CALL timestop(handle)
494 :
495 7240 : END SUBROUTINE geoopt_bfgs
496 :
497 : ! **************************************************************************************************
498 : !> \brief Damped inverse-BFGS update, using only a matvec and local rank-two updates.
499 : !> \param ndf Number of optimization variables
500 : !> \param dx Actual change of optimization variables (not modified by damping)
501 : !> \param dg Change of gradients in the same coordinate system
502 : !> \param hess_mat Inverse Hessian, updated in place
503 : !> \param work Replicated workspace for H*dg
504 : ! **************************************************************************************************
505 502 : SUBROUTINE bfgs_inverse_update(ndf, dx, dg, hess_mat, work)
506 : INTEGER, INTENT(IN) :: ndf
507 : REAL(KIND=dp), INTENT(IN) :: dx(ndf), dg(ndf)
508 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
509 : REAL(KIND=dp), INTENT(OUT) :: work(ndf)
510 :
511 : REAL(KIND=dp), PARAMETER :: min_curvature = 0.2_dp
512 :
513 : INTEGER :: handle, i, ig, invalid, j, jg, nc, nr
514 502 : INTEGER, DIMENSION(:), POINTER :: cols, rows
515 : REAL(KIND=dp) :: coeff, rho, sy, theta, yhy
516 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
517 502 : POINTER :: data
518 502 : REAL(KIND=dp), DIMENSION(ndf) :: secant
519 :
520 502 : CALL timeset("bfgs_inverse_update", handle)
521 157856 : IF (.NOT. ALL(ieee_is_finite(dx)) .OR. .NOT. ALL(ieee_is_finite(dg))) THEN
522 0 : CPABORT("Non-finite displacement or gradient difference in inverse BFGS")
523 : END IF
524 78928 : work = 0.0_dp
525 157354 : IF (NORM2(dx) == 0.0_dp .OR. NORM2(dg) == 0.0_dp) THEN
526 0 : CALL timestop(handle)
527 0 : RETURN
528 : END IF
529 :
530 502 : CALL cp_fm_matvec(hess_mat, dg, work)
531 78928 : sy = DOT_PRODUCT(dx, dg)
532 78928 : yhy = DOT_PRODUCT(dg, work)
533 502 : IF (.NOT. ieee_is_finite(sy) .OR. .NOT. ieee_is_finite(yhy)) THEN
534 0 : CPABORT("Non-finite curvature in inverse BFGS")
535 : END IF
536 502 : IF (yhy <= 0.0_dp) THEN
537 : CALL cp_warn(__LOCATION__, &
538 0 : "BFGS: non-positive inverse curvature; resetting to unity.")
539 0 : CALL cp_fm_set_all(hess_mat, alpha=0.0_dp, beta=1.0_dp)
540 0 : CALL timestop(handle)
541 0 : RETURN
542 : END IF
543 :
544 78928 : secant = dx
545 502 : IF (sy < min_curvature*yhy) THEN
546 : ! Dual-space damping: y^T*s_tilde >= min_curvature*y^T*H*y.
547 : ! Unlike damping y with B*s, this does not need the direct Hessian.
548 4 : theta = ((1.0_dp - min_curvature)*yhy)/(yhy - sy)
549 712 : secant = theta*dx + (1.0_dp - theta)*work
550 712 : sy = DOT_PRODUCT(secant, dg)
551 : END IF
552 502 : IF (sy <= TINY(1.0_dp)/EPSILON(1.0_dp)) THEN
553 0 : CALL timestop(handle)
554 0 : RETURN
555 : END IF
556 502 : rho = 1.0_dp/sy
557 502 : coeff = (1.0_dp + yhy*rho)*rho
558 502 : IF (.NOT. ieee_is_finite(coeff)) THEN
559 : CALL cp_warn(__LOCATION__, &
560 0 : "BFGS: unrepresentable inverse update; retaining the previous matrix.")
561 0 : CALL timestop(handle)
562 0 : RETURN
563 : END IF
564 :
565 : CALL cp_fm_get_info(hess_mat, row_indices=rows, col_indices=cols, &
566 502 : local_data=DATA, nrow_local=nr, ncol_local=nc)
567 502 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(nc,nr,rows,cols,data,coeff,rho,secant,work) PRIVATE(i,j,ig,jg)
568 : DO j = 1, nc
569 : jg = cols(j)
570 : DO i = 1, nr
571 : ig = rows(i)
572 : DATA(i, j) = DATA(i, j) + coeff*(secant(ig)*secant(jg)) - &
573 : rho*(secant(ig)*work(jg) + work(ig)*secant(jg))
574 : END DO
575 : END DO
576 : !$OMP END PARALLEL DO
577 502 : invalid = 0
578 6569539 : IF (.NOT. ALL(ieee_is_finite(DATA(1:nr, 1:nc)))) invalid = 1
579 502 : CALL hess_mat%matrix_struct%para_env%max(invalid)
580 502 : IF (invalid /= 0) THEN
581 : CALL cp_warn(__LOCATION__, &
582 0 : "BFGS: non-finite inverse matrix; resetting to unity.")
583 0 : CALL cp_fm_set_all(hess_mat, alpha=0.0_dp, beta=1.0_dp)
584 : END IF
585 502 : CALL timestop(handle)
586 502 : END SUBROUTINE bfgs_inverse_update
587 :
588 : ! **************************************************************************************************
589 : !> \brief Form -H*g and recover a descent direction if roundoff has corrupted H.
590 : !> \param hess_mat Inverse Hessian
591 : !> \param g Gradient
592 : !> \param dr Search direction
593 : ! **************************************************************************************************
594 512 : SUBROUTINE bfgs_inverse_step(hess_mat, g, dr)
595 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
596 : REAL(KIND=dp), INTENT(IN) :: g(:)
597 : REAL(KIND=dp), INTENT(OUT) :: dr(:)
598 :
599 : REAL(KIND=dp) :: slope
600 :
601 79712 : IF (.NOT. ALL(ieee_is_finite(g))) THEN
602 0 : CALL cp_abort(__LOCATION__, "Non-finite gradient in inverse BFGS")
603 : END IF
604 79712 : dr = 0.0_dp
605 79712 : IF (MAXVAL(ABS(g)) == 0.0_dp) RETURN
606 512 : CALL cp_fm_matvec(hess_mat, g, dr, alpha=-1.0_dp)
607 79712 : slope = DOT_PRODUCT(g, dr)
608 512 : IF (.NOT. ieee_is_finite(slope) .OR. slope >= 0.0_dp) THEN
609 : CALL cp_warn(__LOCATION__, &
610 0 : "BFGS: inverse matrix is not a descent model; resetting to unity.")
611 0 : CALL cp_fm_set_all(hess_mat, alpha=0.0_dp, beta=1.0_dp)
612 512 : dr = -g
613 : END IF
614 : END SUBROUTINE bfgs_inverse_step
615 :
616 : ! **************************************************************************************************
617 : !> \brief Convert a restart matrix once, not during normal inverse-BFGS iterations.
618 : !> \param matrix Input matrix, replaced by its (regularized) inverse
619 : !> \param from_inverse Whether the input already represents the inverse Hessian
620 : ! **************************************************************************************************
621 6 : SUBROUTINE bfgs_invert_matrix(matrix, from_inverse)
622 : TYPE(cp_fm_type), INTENT(IN) :: matrix
623 : LOGICAL, INTENT(IN) :: from_inverse
624 :
625 : INTEGER :: handle, info, ndf
626 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
627 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
628 : TYPE(cp_fm_type) :: eigenvectors, scratch
629 :
630 6 : CALL timeset("bfgs_invert_matrix", handle)
631 6 : CALL cp_fm_get_info(matrix, matrix_struct=matrix_struct, nrow_global=ndf)
632 6 : CALL cp_fm_create(scratch, matrix_struct, name="bfgs_restart_scratch")
633 6 : CALL cp_fm_create(eigenvectors, matrix_struct, name="bfgs_restart_eigenvectors")
634 18 : ALLOCATE (eigenvalues(ndf))
635 6 : CALL cp_fm_to_fm(matrix, scratch)
636 6 : CALL choose_eigv_solver(scratch, eigenvectors, eigenvalues, info=info)
637 6 : IF (info /= 0) THEN
638 0 : CALL cp_abort(__LOCATION__, "BFGS restart matrix diagonalization failed")
639 : END IF
640 246 : IF (.NOT. ALL(ieee_is_finite(eigenvalues))) THEN
641 0 : CALL cp_abort(__LOCATION__, "Non-finite BFGS restart spectrum")
642 : END IF
643 6 : IF (from_inverse) THEN
644 : ! Do not apply the direct-Hessian cutoff to inverse-Hessian eigenvalues.
645 220 : IF (MINVAL(eigenvalues) <= 0.0_dp) THEN
646 : CALL cp_abort(__LOCATION__, &
647 0 : "BFGS inverse restart matrix is not positive definite")
648 : END IF
649 : ELSE
650 26 : IF (MINVAL(eigenvalues) < min_hessian_eigenvalue) THEN
651 : CALL cp_warn(__LOCATION__, &
652 0 : "BFGS: regularizing the direct restart before conversion to INVERSE.")
653 : END IF
654 26 : eigenvalues(:) = MAX(min_hessian_eigenvalue, eigenvalues(:))
655 : END IF
656 246 : eigenvalues(:) = 1.0_dp/eigenvalues(:)
657 246 : IF (.NOT. ALL(ieee_is_finite(eigenvalues))) THEN
658 : CALL cp_abort(__LOCATION__, &
659 0 : "Unrepresentable inverse of BFGS restart matrix")
660 : END IF
661 6 : CALL cp_fm_to_fm(eigenvectors, scratch)
662 6 : CALL cp_fm_column_scale(scratch, eigenvalues)
663 6 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, 1.0_dp, scratch, eigenvectors, 0.0_dp, matrix)
664 6 : CALL cp_fm_release(scratch)
665 6 : CALL cp_fm_release(eigenvectors)
666 6 : DEALLOCATE (eigenvalues)
667 6 : CALL timestop(handle)
668 12 : END SUBROUTINE bfgs_invert_matrix
669 :
670 : ! **************************************************************************************************
671 : !> \brief Read the metadata of a tagged BFGS restart, or identify a legacy raw Hessian.
672 : !> \param unit File unit, returned positioned at the first matrix column
673 : !> \param stored_method Matrix representation in the file
674 : !> \param stored_ndf Matrix dimension stored in a tagged file; zero for legacy files
675 : !> \param stored_type Optimization type stored in a tagged file; zero for legacy files
676 : !> \param tagged Whether the file contains coordinate and representation metadata
677 : !> \param reference_cell Stored CELL_OPT reference cell
678 : !> \param cell_frame Stored canonical-to-input cell frame
679 : ! **************************************************************************************************
680 8 : SUBROUTINE bfgs_read_restart_header(unit, stored_method, stored_ndf, stored_type, tagged, &
681 : reference_cell, cell_frame)
682 : INTEGER, INTENT(IN) :: unit
683 : INTEGER, INTENT(OUT) :: stored_method, stored_ndf, stored_type
684 : LOGICAL, INTENT(OUT) :: tagged
685 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: reference_cell, cell_frame
686 :
687 : CHARACTER(LEN=16) :: magic
688 : INTEGER :: ios, version
689 :
690 8 : stored_method = bfgs_direct
691 8 : stored_ndf = 0
692 8 : stored_type = 0
693 8 : reference_cell = 0.0_dp
694 8 : cell_frame = 0.0_dp
695 8 : magic = ""
696 8 : READ (unit, IOSTAT=ios) magic
697 8 : tagged = ios == 0 .AND. magic == bfgs_restart_magic
698 8 : REWIND (unit)
699 8 : IF (tagged) THEN
700 8 : READ (unit, IOSTAT=ios) magic, version, stored_method, stored_ndf, stored_type, &
701 16 : reference_cell, cell_frame
702 8 : IF (ios /= 0) CPABORT("Invalid BFGS restart header")
703 8 : IF (version /= bfgs_restart_version) THEN
704 0 : CALL cp_abort(__LOCATION__, "Unsupported BFGS restart version")
705 : END IF
706 8 : IF (stored_method /= bfgs_direct .AND. stored_method /= bfgs_inverse) THEN
707 0 : CPABORT("Unknown matrix representation in BFGS restart")
708 : END IF
709 8 : IF (stored_ndf <= 0) CPABORT("Invalid matrix dimension in BFGS restart")
710 208 : IF (.NOT. ALL(ieee_is_finite(reference_cell)) .OR. .NOT. ALL(ieee_is_finite(cell_frame))) THEN
711 0 : CPABORT("Non-finite coordinate frame in BFGS restart")
712 : END IF
713 : END IF
714 8 : END SUBROUTINE bfgs_read_restart_header
715 :
716 : ! **************************************************************************************************
717 : !> \brief Read and validate a tagged or legacy BFGS restart on the I/O rank.
718 : !> \param unit File unit, returned positioned at the first matrix column
719 : !> \param ndf Expected dimension
720 : !> \param opt_type Expected optimization type
721 : !> \param stored_method Matrix representation in the file
722 : !> \param tagged Whether the file contains coordinate and representation metadata
723 : !> \param reference_cell Stored CELL_OPT reference cell
724 : !> \param cell_frame Stored canonical-to-input cell frame
725 : ! **************************************************************************************************
726 4 : SUBROUTINE bfgs_read_header(unit, ndf, opt_type, stored_method, tagged, reference_cell, cell_frame)
727 : INTEGER, INTENT(IN) :: unit, ndf, opt_type
728 : INTEGER, INTENT(OUT) :: stored_method
729 : LOGICAL, INTENT(OUT) :: tagged
730 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: reference_cell, cell_frame
731 :
732 : INTEGER :: ios, j, stored_ndf, stored_type
733 : REAL(KIND=dp) :: extra
734 8 : REAL(KIND=dp), DIMENSION(ndf) :: column
735 :
736 : CALL bfgs_read_restart_header(unit, stored_method, stored_ndf, stored_type, tagged, &
737 4 : reference_cell, cell_frame)
738 4 : IF (tagged) THEN
739 4 : IF (stored_ndf /= ndf) THEN
740 : CALL cp_abort(__LOCATION__, &
741 0 : "BFGS restart dimension does not match the optimization variables")
742 : END IF
743 4 : IF (stored_type /= opt_type) THEN
744 0 : CALL cp_abort(__LOCATION__, "BFGS restart optimization type does not match")
745 : END IF
746 : END IF
747 :
748 : ! Preflight before any collective read: reject truncated/non-finite payloads and
749 : ! incompatible legacy dimensions with a CP2K error, not a rank-local Fortran I/O failure.
750 136 : DO j = 1, ndf
751 132 : READ (unit, IOSTAT=ios) column
752 132 : IF (ios /= 0) THEN
753 : CALL cp_abort(__LOCATION__, &
754 0 : "Truncated BFGS restart or incompatible legacy matrix dimension")
755 : END IF
756 9784 : IF (.NOT. ALL(ieee_is_finite(column))) THEN
757 0 : CALL cp_abort(__LOCATION__, "Non-finite matrix in BFGS restart")
758 : END IF
759 : END DO
760 4 : READ (unit, IOSTAT=ios) extra
761 4 : IF (ios /= iostat_end) THEN
762 : CALL cp_abort(__LOCATION__, &
763 0 : "Extra matrix columns or invalid trailing data in BFGS restart")
764 : END IF
765 4 : REWIND (unit)
766 4 : IF (tagged) READ (unit)
767 4 : END SUBROUTINE bfgs_read_header
768 :
769 : ! **************************************************************************************************
770 : !> \brief Read a common BFGS restart, restoring its coordinates and converting if needed.
771 : !> \param geo_section BFGS input section
772 : !> \param hess_mat Matrix to initialize
773 : !> \param logger Logger
774 : !> \param gopt_env Optimization environment
775 : !> \param x0 Current coordinates, re-expressed in the stored reference cell for CELL_OPT
776 : !> \param matrix_method Requested in-memory representation
777 : ! **************************************************************************************************
778 8 : SUBROUTINE read_bfgs_hessian(geo_section, hess_mat, logger, gopt_env, x0, matrix_method)
779 : TYPE(section_vals_type), POINTER :: geo_section
780 : TYPE(cp_fm_type), INTENT(INOUT) :: hess_mat
781 : TYPE(cp_logger_type), POINTER :: logger
782 : TYPE(gopt_f_type), POINTER :: gopt_env
783 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x0
784 : INTEGER, INTENT(IN) :: matrix_method
785 :
786 : CHARACTER(LEN=default_path_length) :: filename
787 : INTEGER :: ndf, stored_method, unit
788 : LOGICAL :: tagged
789 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_frame, reference_cell
790 : TYPE(mp_para_env_type), POINTER :: para_env
791 :
792 8 : para_env => gopt_env%force_env%para_env
793 8 : ndf = SIZE(x0)
794 8 : unit = -1
795 8 : IF (para_env%is_source()) THEN
796 4 : CALL section_vals_val_get(geo_section, "BFGS%RESTART_FILE_NAME", c_val=filename)
797 4 : IF (LEN_TRIM(filename) == 0) filename = TRIM(logger%iter_info%project_name)//"-BFGS.Hessian"
798 : CALL open_file(file_name=TRIM(filename), file_status="OLD", &
799 4 : file_form="UNFORMATTED", file_action="READ", unit_number=unit)
800 4 : CALL bfgs_read_header(unit, ndf, gopt_env%type_id, stored_method, tagged, reference_cell, cell_frame)
801 : END IF
802 8 : CALL para_env%bcast(stored_method, 0)
803 8 : CALL para_env%bcast(tagged, 0)
804 8 : CALL para_env%bcast(reference_cell, 0)
805 8 : CALL para_env%bcast(cell_frame, 0)
806 8 : CALL cp_fm_read_unformatted(hess_mat, unit)
807 8 : IF (para_env%is_source()) CALL close_file(unit_number=unit)
808 :
809 8 : IF (gopt_env%type_id == default_cell_method_id) THEN
810 2 : IF (tagged) THEN
811 2 : CALL bfgs_restore_cell_frame(gopt_env, x0, reference_cell, cell_frame)
812 : ELSE
813 : CALL cp_warn(__LOCATION__, "Legacy BFGS restart has no reference cell; "// &
814 0 : "assuming the current CELL_OPT coordinate frame. Exact frame recovery is not possible.")
815 : END IF
816 : END IF
817 8 : IF (stored_method /= matrix_method) THEN
818 6 : CALL bfgs_invert_matrix(hess_mat, stored_method == bfgs_inverse)
819 : END IF
820 8 : END SUBROUTINE read_bfgs_hessian
821 :
822 : ! **************************************************************************************************
823 : !> \brief Re-express coordinates in the saved frame without changing physical positions.
824 : !> \param gopt_env Optimization environment
825 : !> \param x0 Optimization coordinates
826 : !> \param reference_cell Reference cell that defines the stored Hessian coordinates
827 : !> \param cell_frame Canonical-to-input cell orientation used when writing the matrix
828 : ! **************************************************************************************************
829 2 : SUBROUTINE bfgs_restore_cell_frame(gopt_env, x0, reference_cell, cell_frame)
830 : TYPE(gopt_f_type), POINTER :: gopt_env
831 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x0
832 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: reference_cell, cell_frame
833 :
834 : INTEGER :: i
835 : REAL(KIND=dp), DIMENSION(3, 3) :: transform
836 :
837 : ! A normal continuation preserves this orientation even when lengths and angles change.
838 : ! Do not silently attach a matrix to a separately rotated or reordered structure.
839 26 : IF (MAXVAL(ABS(cell_frame - gopt_env%cell_env%opt_to_input)) > 1.0E-8_dp) THEN
840 : CALL cp_abort(__LOCATION__, &
841 0 : "BFGS restart cell orientation differs; use a matching structure or RESTART_HESSIAN F")
842 : END IF
843 26 : IF (ABS(det_3x3(reference_cell)) <= EPSILON(1.0_dp)*MAXVAL(ABS(reference_cell))**3) THEN
844 0 : CPABORT("Singular reference cell in BFGS restart")
845 : END IF
846 134 : transform = MATMUL(reference_cell, inv_3x3(gopt_env%h_ref))
847 62 : DO i = 1, SIZE(x0) - 6, 3
848 1142 : x0(i:i + 2) = MATMUL(transform, x0(i:i + 2))
849 : END DO
850 26 : gopt_env%h_ref = reference_cell
851 2 : END SUBROUTINE bfgs_restore_cell_frame
852 :
853 : ! **************************************************************************************************
854 : !> \brief ...
855 : !> \param ndf ...
856 : !> \param dg ...
857 : !> \param eigval ...
858 : !> \param work ...
859 : !> \param eigvec_mat ...
860 : !> \param g ...
861 : !> \param para_env ...
862 : ! **************************************************************************************************
863 1060 : SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
864 :
865 : INTEGER, INTENT(IN) :: ndf
866 : REAL(KIND=dp), INTENT(INOUT) :: dg(ndf), eigval(ndf), work(ndf)
867 : TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat
868 : REAL(KIND=dp), INTENT(INOUT) :: g(ndf)
869 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
870 :
871 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rat_fun_opt'
872 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
873 :
874 : INTEGER :: handle, i, indf, iref, iter, j, k, l, &
875 : maxit, ncol_local, nrow_local
876 530 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
877 : LOGICAL :: bisec, conv, fail, set
878 : REAL(KIND=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
879 : ln, lp, ssize, step, stol
880 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
881 530 : POINTER :: local_data
882 :
883 530 : CALL timeset(routineN, handle)
884 :
885 530 : stol = 1.0E-8_dp
886 530 : ssize = 0.2_dp
887 530 : maxit = 999
888 530 : fail = .FALSE.
889 530 : bisec = .FALSE.
890 :
891 70634 : dg = 0._dp
892 :
893 : CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
894 530 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
895 :
896 40652 : DO i = 1, nrow_local
897 40122 : j = row_indices(i)
898 10898936 : DO k = 1, ncol_local
899 10858284 : l = col_indices(k)
900 10898406 : dg(l) = dg(l) + local_data(i, k)*g(j)
901 : END DO
902 : END DO
903 530 : CALL para_env%sum(dg)
904 :
905 530 : set = .FALSE.
906 :
907 : DO
908 :
909 : ! calculating Lambda
910 :
911 554 : lp = 0.0_dp
912 554 : iref = 1
913 554 : ln = 0.0_dp
914 554 : IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
915 :
916 554 : conv = .FALSE.
917 554 : iter = 0
918 : DO
919 1260 : iter = iter + 1
920 1260 : fun = 0.0_dp
921 1260 : fung = 0.0_dp
922 153612 : DO indf = 1, ndf
923 152352 : fun = fun + dg(indf)**2/(ln - eigval(indf))
924 153612 : fung = fung - dg(indf)**2/((ln - eigval(indf))**2)
925 : END DO
926 1260 : fun = fun - ln
927 1260 : fung = fung - one
928 1260 : step = fun/fung
929 1260 : ln = ln - step
930 1260 : IF (ABS(step) < stol) THEN
931 : conv = .TRUE.
932 : EXIT
933 : END IF
934 706 : IF (iter >= maxit) EXIT
935 : END DO
936 : outer: DO
937 578 : IF (.NOT. conv) THEN
938 744 : bisec = .TRUE.
939 744 : iter = 0
940 744 : maxit = 9999
941 744 : lam1 = 0.0_dp
942 744 : IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
943 : fun1 = 0.0_dp
944 744 : DO indf = 1, ndf
945 744 : fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
946 : END DO
947 24 : fun1 = fun1 - lam1
948 24 : step = ABS(lam1)/1000.0_dp
949 : IF (step < ssize) step = ssize
950 : inner: DO
951 240000 : iter = iter + 1
952 240000 : IF (iter > maxit) THEN
953 : ln = 0.0_dp
954 554 : lp = 0.0_dp
955 : fail = .TRUE.
956 : EXIT outer
957 : END IF
958 239976 : fun2 = 0.0_dp
959 239976 : lam2 = lam1 - iter*step
960 7439256 : DO indf = 1, ndf
961 7439256 : fun2 = fun2 + dg(indf)**2/(lam2 - eigval(indf))
962 : END DO
963 239976 : fun2 = fun2 - lam2
964 239976 : IF (fun2*fun1 < 0.0_dp) THEN
965 : iter = 0
966 : DO
967 0 : iter = iter + 1
968 0 : IF (iter > maxit) THEN
969 : ln = 0.0_dp
970 : lp = 0.0_dp
971 : fail = .TRUE.
972 : EXIT outer
973 : END IF
974 0 : step = (lam1 + lam2)/2
975 0 : fun3 = 0.0_dp
976 0 : DO indf = 1, ndf
977 0 : fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
978 : END DO
979 0 : fun3 = fun3 - step
980 :
981 0 : IF (ABS(step - lam2) < stol) THEN
982 : ln = step
983 : EXIT inner
984 : END IF
985 :
986 0 : IF (fun3*fun1 < stol) THEN
987 : lam2 = step
988 : ELSE
989 0 : lam1 = step
990 : END IF
991 : END DO
992 : END IF
993 : END DO inner
994 : END IF
995 554 : IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
996 : (eigval(iref) > 0.0_dp))) THEN
997 :
998 24 : IF (.NOT. bisec) THEN
999 : conv = .FALSE.
1000 : CYCLE outer
1001 : END IF
1002 : ln = 0.0_dp
1003 : lp = 0.0_dp
1004 : fail = .TRUE.
1005 : END IF
1006 : EXIT outer
1007 : END DO outer
1008 :
1009 554 : IF (fail .AND. .NOT. set) THEN
1010 744 : set = .TRUE.
1011 744 : DO indf = 1, ndf
1012 744 : eigval(indf) = eigval(indf)*work(indf)
1013 : END DO
1014 : CYCLE
1015 : END IF
1016 :
1017 530 : IF (.NOT. set) THEN
1018 69890 : work(1:ndf) = one
1019 : END IF
1020 :
1021 70634 : DO indf = 1, ndf
1022 70634 : eigval(indf) = eigval(indf) - ln
1023 : END DO
1024 : EXIT
1025 : END DO
1026 :
1027 530 : CALL timestop(handle)
1028 :
1029 530 : END SUBROUTINE rat_fun_opt
1030 :
1031 : ! **************************************************************************************************
1032 : !> \brief Update the direct Hessian with the damped BFGS formula
1033 : !> \param ndf Number of optimization variables
1034 : !> \param dx Change of the optimization variables
1035 : !> \param dg Change of the gradients in the same coordinate system
1036 : !> \param hess_mat Direct Hessian, updated in place
1037 : !> \param work Replicated workspace, holds B*dx on exit
1038 : !> \param para_env ...
1039 : !> \param output_unit Geometry log unit, activates at the default MEDIUM print level
1040 : !> \param hess_reset Set if the update had to be refused and the caller must
1041 : !> collectively rebuild the initial Hessian
1042 : ! **************************************************************************************************
1043 5970 : SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env, output_unit, hess_reset)
1044 : INTEGER, INTENT(IN) :: ndf
1045 : REAL(KIND=dp), INTENT(INOUT) :: dx(ndf), dg(ndf)
1046 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1047 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
1048 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
1049 : INTEGER, INTENT(IN) :: output_unit
1050 : LOGICAL, INTENT(OUT) :: hess_reset
1051 :
1052 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bfgs'
1053 : INTEGER, PARAMETER :: exit_reset = 2, exit_skip = 1, &
1054 : exit_update = 0
1055 : REAL(KIND=dp), PARAMETER :: min_curvature = 0.2_dp, one = 1.0_dp, &
1056 : min_sbs = 10.0_dp*TINY(one), &
1057 : zero = 0.0_dp
1058 :
1059 : INTEGER :: exit_flag, handle, i, invalid, j, k, l, &
1060 : ncol_local, nrow_local
1061 2985 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1062 : REAL(KIND=dp) :: dxw, eta_i, eta_j, gdx, sbs, seta, sy, &
1063 : theta
1064 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1065 2985 : POINTER :: local_hes
1066 :
1067 2985 : CALL timeset(routineN, handle)
1068 2985 : hess_reset = .FALSE.
1069 :
1070 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1071 2985 : local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
1072 :
1073 : ! Refuse corrupted input instead of poisoning the Hessian with Inf/NaN.
1074 494574 : IF (.NOT. ALL(ieee_is_finite(dx(1:ndf))) .OR. .NOT. ALL(ieee_is_finite(dg(1:ndf)))) THEN
1075 0 : CPABORT("Non-finite displacement or gradient difference")
1076 : END IF
1077 :
1078 247287 : work = zero
1079 145209 : DO i = 1, nrow_local
1080 142224 : j = row_indices(i)
1081 23215557 : DO k = 1, ncol_local
1082 23070348 : l = col_indices(k)
1083 23212572 : work(j) = work(j) + local_hes(i, k)*dx(l)
1084 : END DO
1085 : END DO
1086 :
1087 2985 : IF (PRESENT(para_env)) CALL para_env%sum(work)
1088 :
1089 247287 : sy = DOT_PRODUCT(dg(1:ndf), dx(1:ndf))
1090 247287 : sbs = DOT_PRODUCT(dx(1:ndf), work(1:ndf))
1091 :
1092 : ! One collective decision covers every early exit, so that all ranks
1093 : ! take the same path. sbs and sy share the units of the energy, so the
1094 : ! curvature floor is relative to sy. min_sbs only fences off values for
1095 : ! which the divisions in the update below would overflow: with a
1096 : ! smaller sbs, 0.2*sbs (the least possible seta) falls below the
1097 : ! smallest normal number. The floor is never negative, so a non-positive
1098 : ! sbs (lost positive definiteness) is caught as well.
1099 2985 : exit_flag = exit_update
1100 494574 : IF (DOT_PRODUCT(dx(1:ndf), dx(1:ndf)) == zero .OR. &
1101 : DOT_PRODUCT(dg(1:ndf), dg(1:ndf)) == zero) THEN
1102 0 : exit_flag = exit_skip
1103 2985 : ELSE IF (.NOT. ieee_is_finite(sy) .OR. .NOT. ieee_is_finite(sbs) .OR. &
1104 : sbs <= MAX(EPSILON(one)*ABS(sy), min_sbs)) THEN
1105 0 : exit_flag = exit_reset
1106 : END IF
1107 : ! The max folds the votes into the most conservative action all ranks
1108 : ! can agree on. The exit codes are ordered by severity, so a skip vote
1109 : ! on one rank and a reset vote on another reset both. Every rank takes
1110 : ! the same branch, so the rebuild stays collective.
1111 2985 : IF (PRESENT(para_env)) CALL para_env%max(exit_flag)
1112 :
1113 2985 : IF (exit_flag == exit_skip) THEN
1114 : ! A zero displacement or gradient difference carries no curvature
1115 : ! information, so there is nothing to learn from this step.
1116 2985 : ELSE IF (exit_flag == exit_reset) THEN
1117 : ! The Hessian has lost positive definiteness along dx or the predicted
1118 : ! curvature is not representable, so the secant pair cannot be damped:
1119 : ! the caller restarts from the initial Hessian and skips this update.
1120 : CALL cp_warn(__LOCATION__, &
1121 0 : "BFGS: vanishing or non-positive predicted curvature, resetting the Hessian")
1122 0 : hess_reset = .TRUE.
1123 : ELSE
1124 : ! Damped update after M.J.D. Powell, Lecture Notes in Mathematics 630 (1978),
1125 : ! eq. 3.7: whenever the measured curvature sy drops below 20% of the predicted
1126 : ! curvature s.B.s, replace y by eta = theta*y + (1-theta)*B*s. The blended
1127 : ! curvature satisfies eta.s = 0.2*s.B.s exactly, so the update stays positive
1128 : ! definite in exact arithmetic, also for secant pairs that no line search
1129 : ! checked. Because sbs passed the floor above, seta >= 0.2*sbs stays a
1130 : ! normal number, so one/seta cannot overflow and no further floor on seta
1131 : ! is needed.
1132 2985 : theta = one
1133 2985 : IF (sy < min_curvature*sbs) THEN
1134 65 : theta = ((one - min_curvature)*sbs)/(sbs - sy)
1135 65 : IF (output_unit > 0) THEN
1136 : WRITE (UNIT=output_unit, FMT="(T2,A,F10.4,A,F3.1,A)") &
1137 26 : "BFGS| Damped update, curvature ratio y.s/(s.B.s): ", sy/sbs, &
1138 52 : " (threshold ", min_curvature, ")"
1139 26 : CALL m_flush(output_unit)
1140 : END IF
1141 : END IF
1142 2985 : seta = theta*sy + (one - theta)*sbs
1143 2985 : gdx = one/seta
1144 2985 : dxw = one/sbs
1145 :
1146 145209 : DO i = 1, nrow_local
1147 142224 : j = row_indices(i)
1148 142224 : eta_i = theta*dg(j) + (one - theta)*work(j)
1149 23215557 : DO k = 1, ncol_local
1150 23070348 : l = col_indices(k)
1151 23070348 : eta_j = theta*dg(l) + (one - theta)*work(l)
1152 : local_hes(i, k) = local_hes(i, k) + gdx*eta_i*eta_j - &
1153 23212572 : dxw*work(j)*work(l)
1154 : END DO
1155 : END DO
1156 :
1157 : ! Collective verdict: if any local block is corrupted, every rank resets.
1158 2985 : invalid = 0
1159 23317635 : IF (.NOT. ALL(ieee_is_finite(local_hes(1:nrow_local, 1:ncol_local)))) invalid = 1
1160 2985 : IF (PRESENT(para_env)) CALL para_env%max(invalid)
1161 2985 : IF (invalid /= 0) THEN
1162 0 : CALL cp_warn(__LOCATION__, "BFGS: non-finite Hessian after the update, resetting")
1163 0 : hess_reset = .TRUE.
1164 : END IF
1165 : END IF
1166 :
1167 2985 : CALL timestop(handle)
1168 :
1169 2985 : END SUBROUTINE bfgs
1170 :
1171 : ! **************************************************************************************************
1172 : !> \brief Build the initial direct or inverse Hessian: sanitized model Hessian or unity
1173 : !> \param hess_mat BFGS matrix (direct or inverse Hessian), replaced in place
1174 : !> \param hess_tmp Eigensolver scratch matrix, same structure as hess_mat
1175 : !> \param eigvec Eigensolver eigenvector matrix, same structure as hess_mat
1176 : !> \param eigval Eigensolver workspace of length ndf
1177 : !> \param force_env Environment providing the model Hessian at the current geometry
1178 : !> \param use_model Build the sanitized model Hessian; otherwise fall back to unity
1179 : !> \param inverse Store the inverse of the sanitized model Hessian
1180 : !> \note Must be called by all ranks together: construct_initial_hess reduces internally.
1181 : ! **************************************************************************************************
1182 901 : SUBROUTINE build_initial_hess(hess_mat, hess_tmp, eigvec, eigval, force_env, use_model, inverse)
1183 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat, hess_tmp, eigvec
1184 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigval
1185 : TYPE(force_env_type), POINTER :: force_env
1186 : LOGICAL, INTENT(IN) :: use_model, inverse
1187 :
1188 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_initial_hess'
1189 : REAL(KIND=dp), PARAMETER :: min_eig = 0.1_dp, one = 1.0_dp, &
1190 : zero = 0.0_dp
1191 :
1192 : INTEGER :: handle, indf, info, ndf
1193 :
1194 901 : CALL timeset(routineN, handle)
1195 :
1196 901 : IF (use_model) THEN
1197 767 : CPASSERT(ASSOCIATED(force_env))
1198 767 : ndf = SIZE(eigval)
1199 : ! The raw model Hessian is not guaranteed positive definite, so
1200 : ! eigenfactorize, floor the eigenvalues, and rebuild.
1201 767 : CALL cp_fm_set_all(hess_mat, alpha=zero)
1202 767 : CALL construct_initial_hess(force_env, hess_mat)
1203 767 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
1204 767 : CALL choose_eigv_solver(hess_tmp, eigvec, eigval, info=info)
1205 767 : IF (info /= 0) THEN
1206 : ! In rare cases the diagonalization fails (bug in scalapack?)
1207 0 : CALL cp_warn(__LOCATION__, "BFGS: model Hessian diagonalization failed, using unity")
1208 0 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
1209 0 : eigval = one
1210 : ELSE
1211 16715 : DO indf = 1, ndf
1212 15948 : IF (eigval(indf) < min_eig) eigval(indf) = min_eig
1213 16715 : IF (inverse) eigval(indf) = one/eigval(indf)
1214 : END DO
1215 767 : CALL cp_fm_to_fm(eigvec, hess_tmp)
1216 767 : CALL cp_fm_column_scale(eigvec, eigval)
1217 767 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec, zero, hess_mat)
1218 : END IF
1219 : ELSE
1220 : ! No sanitized prior exists (user opt-out, cell optimization, or
1221 : ! core-shell model): fall back to unity. eigval mirrors the unity
1222 : ! spectrum, keeping the INTENT(OUT) argument defined on every path.
1223 134 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
1224 8558 : eigval = one
1225 : END IF
1226 :
1227 901 : CALL timestop(handle)
1228 :
1229 901 : END SUBROUTINE build_initial_hess
1230 :
1231 : ! **************************************************************************************************
1232 : !> \brief ...
1233 : !> \param ndf ...
1234 : !> \param eigval ...
1235 : !> \param work ...
1236 : ! **************************************************************************************************
1237 530 : SUBROUTINE set_hes_eig(ndf, eigval, work)
1238 : INTEGER, INTENT(IN) :: ndf
1239 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf), work(ndf)
1240 :
1241 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_hes_eig'
1242 : REAL(KIND=dp), PARAMETER :: max_eig = 1000.0_dp, max_neg = -0.5_dp, &
1243 : max_pos = 5.0_dp, min_eig = 0.005_dp, &
1244 : one = 1.0_dp
1245 :
1246 : INTEGER :: handle, indf
1247 :
1248 530 : CALL timeset(routineN, handle)
1249 :
1250 70634 : DO indf = 1, ndf
1251 70634 : IF (eigval(indf) > max_eig) eigval(indf) = max_eig
1252 : END DO
1253 70634 : DO indf = 1, ndf
1254 70634 : IF (eigval(indf) < 0.0_dp) THEN
1255 24 : IF (eigval(indf) < max_neg) THEN
1256 0 : eigval(indf) = max_neg
1257 24 : ELSE IF (eigval(indf) > -min_eig) THEN
1258 24 : eigval(indf) = -min_eig
1259 : END IF
1260 70080 : ELSE IF (eigval(indf) < 1000.0_dp) THEN
1261 70056 : IF (eigval(indf) < min_eig) THEN
1262 168 : eigval(indf) = min_eig
1263 69888 : ELSE IF (eigval(indf) > max_pos) THEN
1264 5 : eigval(indf) = max_pos
1265 : END IF
1266 : END IF
1267 : END DO
1268 :
1269 70634 : DO indf = 1, ndf
1270 70634 : IF (eigval(indf) < 0.0_dp) THEN
1271 24 : work(indf) = -one
1272 : ELSE
1273 70080 : work(indf) = one
1274 : END IF
1275 : END DO
1276 :
1277 530 : CALL timestop(handle)
1278 :
1279 530 : END SUBROUTINE set_hes_eig
1280 :
1281 : ! **************************************************************************************************
1282 : !> \brief ...
1283 : !> \param ndf ...
1284 : !> \param eigval ...
1285 : !> \param eigvec_mat ...
1286 : !> \param hess_tmp ...
1287 : !> \param dr ...
1288 : !> \param g ...
1289 : !> \param para_env ...
1290 : !> \param use_rfo ...
1291 : ! **************************************************************************************************
1292 11610 : SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
1293 :
1294 : INTEGER, INTENT(IN) :: ndf
1295 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf)
1296 : TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat, hess_tmp
1297 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
1298 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
1299 : LOGICAL :: use_rfo
1300 :
1301 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1302 :
1303 : INTEGER :: i, indf, j, k, l, ncol_local, nrow_local
1304 3870 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1305 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1306 3870 : POINTER :: local_data
1307 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1308 : TYPE(cp_fm_type) :: tmp
1309 :
1310 3870 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
1311 3870 : IF (use_rfo) THEN
1312 70634 : DO indf = 1, ndf
1313 70634 : eigval(indf) = one/eigval(indf)
1314 : END DO
1315 : ELSE
1316 201208 : DO indf = 1, ndf
1317 201208 : eigval(indf) = one/MAX(min_hessian_eigenvalue, eigval(indf))
1318 : END DO
1319 : END IF
1320 :
1321 3870 : CALL cp_fm_column_scale(hess_tmp, eigval)
1322 3870 : CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
1323 3870 : CALL cp_fm_create(tmp, matrix_struct, name="tmp")
1324 3870 : CALL cp_fm_set_all(tmp, alpha=zero)
1325 :
1326 3870 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
1327 :
1328 3870 : CALL cp_fm_transpose(tmp, hess_tmp)
1329 3870 : CALL cp_fm_release(tmp)
1330 :
1331 : ! New step
1332 :
1333 : CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
1334 3870 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1335 :
1336 271842 : dr = 0.0_dp
1337 159654 : DO i = 1, nrow_local
1338 155784 : j = row_indices(i)
1339 25019346 : DO k = 1, ncol_local
1340 24859692 : l = col_indices(k)
1341 25015476 : dr(j) = dr(j) - local_data(i, k)*g(l)
1342 : END DO
1343 : END DO
1344 :
1345 3870 : CALL para_env%sum(dr)
1346 :
1347 3870 : END SUBROUTINE geoopt_get_step
1348 :
1349 : ! **************************************************************************************************
1350 : !> \brief ...
1351 : !> \param ndf ...
1352 : !> \param step ...
1353 : !> \param rad ...
1354 : !> \param rat ...
1355 : !> \param dr ...
1356 : !> \param output_unit ...
1357 : ! **************************************************************************************************
1358 3870 : SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
1359 : INTEGER, INTENT(IN) :: ndf
1360 : REAL(KIND=dp), INTENT(INOUT) :: step, rad, rat, dr(ndf)
1361 : INTEGER, INTENT(IN) :: output_unit
1362 :
1363 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trust_radius'
1364 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
1365 :
1366 : INTEGER :: handle
1367 : REAL(KIND=dp) :: scal
1368 :
1369 3870 : CALL timeset(routineN, handle)
1370 :
1371 271842 : step = MAXVAL(ABS(dr))
1372 3870 : scal = MAX(one, rad/step)
1373 :
1374 3870 : IF (step > rad) THEN
1375 370 : rat = rad/step
1376 370 : CALL DSCAL(ndf, rat, dr, 1)
1377 370 : step = rad
1378 370 : IF (output_unit > 0) THEN
1379 : WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") &
1380 195 : " Step is scaled; Scaling factor = ", rat
1381 195 : CALL m_flush(output_unit)
1382 : END IF
1383 : END IF
1384 3870 : CALL timestop(handle)
1385 :
1386 3870 : END SUBROUTINE trust_radius
1387 :
1388 : ! **************************************************************************************************
1389 : !> \brief ...
1390 : !> \param ndf ...
1391 : !> \param work ...
1392 : !> \param hess_mat ...
1393 : !> \param dr ...
1394 : !> \param g ...
1395 : !> \param pred ...
1396 : !> \param para_env ...
1397 : ! **************************************************************************************************
1398 7740 : SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, pred, para_env)
1399 :
1400 : INTEGER, INTENT(IN) :: ndf
1401 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
1402 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1403 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf), pred
1404 : TYPE(mp_para_env_type), POINTER :: para_env
1405 :
1406 : CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_predict'
1407 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
1408 :
1409 : INTEGER :: handle, i, j, k, l, ncol_local, &
1410 : nrow_local
1411 3870 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1412 : REAL(KIND=dp) :: ener1, ener2
1413 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1414 3870 : POINTER :: local_data
1415 :
1416 3870 : CALL timeset(routineN, handle)
1417 :
1418 271842 : ener1 = DOT_PRODUCT(g(1:ndf), dr(1:ndf))
1419 :
1420 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1421 3870 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1422 :
1423 271842 : work = zero
1424 159654 : DO i = 1, nrow_local
1425 155784 : j = row_indices(i)
1426 25019346 : DO k = 1, ncol_local
1427 24859692 : l = col_indices(k)
1428 25015476 : work(j) = work(j) + local_data(i, k)*dr(l)
1429 : END DO
1430 : END DO
1431 :
1432 3870 : CALL para_env%sum(work)
1433 271842 : ener2 = DOT_PRODUCT(dr(1:ndf), work(1:ndf))
1434 3870 : pred = ener1 + 0.5_dp*ener2
1435 3870 : CALL timestop(handle)
1436 :
1437 3870 : END SUBROUTINE energy_predict
1438 :
1439 : ! **************************************************************************************************
1440 : !> \brief ...
1441 : !> \param rat ...
1442 : !> \param rad ...
1443 : !> \param step ...
1444 : !> \param ediff ...
1445 : ! **************************************************************************************************
1446 505 : SUBROUTINE update_trust_rad(rat, rad, step, ediff)
1447 :
1448 : REAL(KIND=dp), INTENT(INOUT) :: rat, rad, step, ediff
1449 :
1450 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_trust_rad'
1451 : REAL(KIND=dp), PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp
1452 :
1453 : INTEGER :: handle
1454 :
1455 505 : CALL timeset(routineN, handle)
1456 :
1457 505 : IF (rat > 4.0_dp) THEN
1458 0 : IF (ediff < 0.0_dp) THEN
1459 0 : rad = step*0.5_dp
1460 : ELSE
1461 0 : rad = step*0.25_dp
1462 : END IF
1463 505 : ELSE IF (rat > 2.0_dp) THEN
1464 0 : IF (ediff < 0.0_dp) THEN
1465 0 : rad = step*0.75_dp
1466 : ELSE
1467 0 : rad = step*0.5_dp
1468 : END IF
1469 505 : ELSE IF (rat > 4.0_dp/3.0_dp) THEN
1470 0 : IF (ediff < 0.0_dp) THEN
1471 0 : rad = step
1472 : ELSE
1473 0 : rad = step*0.75_dp
1474 : END IF
1475 505 : ELSE IF (rat > 10.0_dp/9.0_dp) THEN
1476 0 : IF (ediff < 0.0_dp) THEN
1477 0 : rad = step*1.25_dp
1478 : ELSE
1479 0 : rad = step
1480 : END IF
1481 505 : ELSE IF (rat > 0.9_dp) THEN
1482 41 : IF (ediff < 0.0_dp) THEN
1483 41 : rad = step*1.5_dp
1484 : ELSE
1485 0 : rad = step*1.25_dp
1486 : END IF
1487 464 : ELSE IF (rat > 0.75_dp) THEN
1488 114 : IF (ediff < 0.0_dp) THEN
1489 109 : rad = step*1.25_dp
1490 : ELSE
1491 5 : rad = step
1492 : END IF
1493 350 : ELSE IF (rat > 0.5_dp) THEN
1494 42 : IF (ediff < 0.0_dp) THEN
1495 42 : rad = step
1496 : ELSE
1497 0 : rad = step*0.75_dp
1498 : END IF
1499 308 : ELSE IF (rat > 0.25_dp) THEN
1500 4 : IF (ediff < 0.0_dp) THEN
1501 4 : rad = step*0.75_dp
1502 : ELSE
1503 0 : rad = step*0.5_dp
1504 : END IF
1505 304 : ELSE IF (ediff < 0.0_dp) THEN
1506 302 : rad = step*0.5_dp
1507 : ELSE
1508 2 : rad = step*0.25_dp
1509 : END IF
1510 :
1511 505 : rad = MAX(rad, min_trust)
1512 505 : rad = MIN(rad, max_trust)
1513 505 : CALL timestop(handle)
1514 :
1515 505 : END SUBROUTINE update_trust_rad
1516 :
1517 : ! **************************************************************************************************
1518 :
1519 : ! **************************************************************************************************
1520 : !> \brief ...
1521 : !> \param geo_section ...
1522 : !> \param hess_mat ...
1523 : !> \param logger ...
1524 : !> \param gopt_env Optimization coordinate system
1525 : !> \param matrix_method Stored matrix representation
1526 : ! **************************************************************************************************
1527 8134 : SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger, gopt_env, matrix_method)
1528 :
1529 : TYPE(section_vals_type), POINTER :: geo_section
1530 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1531 : TYPE(cp_logger_type), POINTER :: logger
1532 : TYPE(gopt_f_type), POINTER :: gopt_env
1533 : INTEGER, INTENT(IN) :: matrix_method
1534 :
1535 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bfgs_hessian'
1536 :
1537 : INTEGER :: handle, hesunit, ndf
1538 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_frame
1539 :
1540 4067 : CALL timeset(routineN, handle)
1541 4067 : CALL cp_fm_get_info(hess_mat, nrow_global=ndf)
1542 4067 : cell_frame = 0.0_dp
1543 11435 : IF (gopt_env%type_id == default_cell_method_id) cell_frame = gopt_env%cell_env%opt_to_input
1544 : hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", &
1545 : extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", &
1546 4067 : file_position="REWIND")
1547 4067 : IF (hesunit > 0) THEN
1548 2355 : WRITE (hesunit) bfgs_restart_magic, bfgs_restart_version, matrix_method, ndf, &
1549 32970 : gopt_env%type_id, gopt_env%h_ref, cell_frame
1550 : END IF
1551 : ! Store the native representation. Conversions occur only when reading a different METHOD.
1552 4067 : CALL cp_fm_write_unformatted(hess_mat, hesunit)
1553 4067 : CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
1554 4067 : CALL timestop(handle)
1555 :
1556 4067 : END SUBROUTINE write_bfgs_hessian
1557 :
1558 : ! **************************************************************************************************
1559 : !> \brief Constructs model Hessian as described in https://doi.org/10.1016/0009-2614(95)00646-L.
1560 : !> \param force_env ...
1561 : !> \param hess_mat ...
1562 : !> \author Florian Schiffmann
1563 : ! **************************************************************************************************
1564 767 : SUBROUTINE construct_initial_hess(force_env, hess_mat)
1565 :
1566 : TYPE(force_env_type), POINTER :: force_env
1567 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1568 :
1569 : INTEGER :: i, iat_col, iat_row, iglobal, iind, j, &
1570 : jat_row, jglobal, jind, k, natom, &
1571 : ncol_local, nrow_local, z
1572 767 : INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row
1573 767 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1574 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij
1575 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij
1576 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1577 767 : POINTER :: fixed, local_data
1578 : REAL(KIND=dp), DIMENSION(3, 3) :: alpha, r0
1579 : TYPE(cell_type), POINTER :: cell
1580 : TYPE(cp_subsys_type), POINTER :: subsys
1581 : TYPE(particle_list_type), POINTER :: particles
1582 :
1583 767 : CALL cite_reference(Lindh1995)
1584 :
1585 767 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1586 : CALL cp_subsys_get(subsys, &
1587 767 : particles=particles)
1588 :
1589 767 : alpha = lindh_alpha
1590 767 : r0 = lindh_r0
1591 :
1592 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1593 767 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1594 767 : natom = particles%n_els
1595 2301 : ALLOCATE (at_row(natom))
1596 3068 : ALLOCATE (rho_ij(natom, natom))
1597 2301 : ALLOCATE (d_ij(natom, natom))
1598 3835 : ALLOCATE (r_ij(natom, natom, 3))
1599 2301 : ALLOCATE (fixed(3, natom))
1600 22031 : fixed = 1.0_dp
1601 767 : CALL fix_atom_control(force_env, fixed)
1602 3068 : DO i = 1, 3
1603 34964 : CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1604 : END DO
1605 767 : rho_ij = 0
1606 6083 : at_row = 3
1607 6083 : DO i = 1, natom
1608 5316 : CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
1609 5316 : IF (z <= 10) at_row(i) = 2
1610 11399 : IF (z <= 2) at_row(i) = 1
1611 : END DO
1612 5316 : DO i = 2, natom
1613 4549 : iat_row = at_row(i)
1614 68083 : DO j = 1, i - 1
1615 62767 : jat_row = at_row(j)
1616 : !pbc for a distance vector
1617 251068 : r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
1618 251068 : r_ij(i, j, :) = -r_ij(j, i, :)
1619 251068 : d_ij(j, i) = NORM2(r_ij(j, i, :))
1620 62767 : d_ij(i, j) = d_ij(j, i)
1621 62767 : rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1622 67316 : rho_ij(i, j) = rho_ij(j, i)
1623 : END DO
1624 : END DO
1625 16715 : DO i = 1, ncol_local
1626 15948 : iglobal = col_indices(i)
1627 15948 : iind = MOD(iglobal - 1, 3) + 1
1628 15948 : iat_col = (iglobal + 2)/3
1629 15948 : IF (iat_col > natom) CYCLE
1630 689420 : DO j = 1, nrow_local
1631 672705 : jglobal = row_indices(j)
1632 672705 : jind = MOD(jglobal - 1, 3) + 1
1633 672705 : iat_row = (jglobal + 2)/3
1634 672705 : IF (iat_row > natom) CYCLE
1635 672705 : IF (iat_row /= iat_col) THEN
1636 643608 : IF (d_ij(iat_row, iat_col) < lindh_cutoff) THEN
1637 : local_data(j, i) = local_data(j, i) + &
1638 213120 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1639 : END IF
1640 : ELSE
1641 : local_data(j, i) = local_data(j, i) + &
1642 29097 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1643 : END IF
1644 672705 : IF (iat_col /= iat_row) THEN
1645 643608 : IF (d_ij(iat_row, iat_col) < lindh_cutoff) THEN
1646 : local_data(j, i) = local_data(j, i) - &
1647 : dist_second_deriv(r_ij(iat_col, iat_row, :), &
1648 852480 : iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1649 : END IF
1650 : ELSE
1651 701802 : DO k = 1, natom
1652 672705 : IF (k == iat_col) CYCLE
1653 672705 : IF (d_ij(iat_row, k) < lindh_cutoff) THEN
1654 : local_data(j, i) = local_data(j, i) + &
1655 : dist_second_deriv(r_ij(iat_col, k, :), &
1656 852480 : iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1657 : END IF
1658 : END DO
1659 : END IF
1660 688653 : IF (fixed(jind, iat_row) < 0.5_dp .OR. fixed(iind, iat_col) < 0.5_dp) THEN
1661 10161 : local_data(j, i) = 0.0_dp
1662 10161 : IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
1663 : END IF
1664 : END DO
1665 : END DO
1666 767 : DEALLOCATE (fixed)
1667 767 : DEALLOCATE (rho_ij)
1668 767 : DEALLOCATE (d_ij)
1669 767 : DEALLOCATE (r_ij)
1670 767 : DEALLOCATE (at_row)
1671 :
1672 1534 : END SUBROUTINE construct_initial_hess
1673 :
1674 : ! **************************************************************************************************
1675 : !> \brief ...
1676 : !> \param r1 ...
1677 : !> \param i ...
1678 : !> \param j ...
1679 : !> \param d ...
1680 : !> \param rho ...
1681 : !> \return ...
1682 : ! **************************************************************************************************
1683 426240 : PURE FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv)
1684 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r1
1685 : INTEGER, INTENT(IN) :: i, j
1686 : REAL(KIND=dp), INTENT(IN) :: d, rho
1687 : REAL(KIND=dp) :: deriv
1688 :
1689 426240 : deriv = lindh_stretch_factor*rho*(r1(i)*r1(j))/d**2
1690 426240 : END FUNCTION dist_second_deriv
1691 :
1692 : ! **************************************************************************************************
1693 : !> \brief ...
1694 : !> \param r_ij ...
1695 : !> \param d_ij ...
1696 : !> \param rho_ij ...
1697 : !> \param idir ...
1698 : !> \param jdir ...
1699 : !> \param iat_der ...
1700 : !> \param jat_der ...
1701 : !> \param natom ...
1702 : !> \return ...
1703 : ! **************************************************************************************************
1704 242217 : FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv)
1705 : REAL(KIND=dp), DIMENSION(:, :, :) :: r_ij
1706 : REAL(KIND=dp), DIMENSION(:, :) :: d_ij, rho_ij
1707 : INTEGER :: idir, jdir, iat_der, jat_der, natom
1708 : REAL(KIND=dp) :: deriv
1709 :
1710 : REAL(KIND=dp), PARAMETER :: rho_thre = 1.0E-5_dp
1711 :
1712 : INTEGER :: i, iat, idr, j, jat, jdr
1713 : REAL(KIND=dp) :: d12, d23, d31, D_mat(3, 2), denom1, &
1714 : denom2, denom3, ka1, ka2, ka3, rho12, &
1715 : rho23, rho31, rsst1, rsst2, rsst3
1716 : REAL(KIND=dp), DIMENSION(3) :: r12, r23, r31
1717 :
1718 242217 : deriv = 0._dp
1719 242217 : IF (iat_der == jat_der) THEN
1720 672705 : DO i = 1, natom - 1
1721 643608 : IF (rho_ij(iat_der, i) < rho_thre) CYCLE
1722 3569778 : DO j = i + 1, natom
1723 3344301 : IF (rho_ij(iat_der, j) < rho_thre) CYCLE
1724 968058 : IF (i == iat_der .OR. j == iat_der) CYCLE
1725 968058 : IF (iat_der < i .OR. iat_der > j) THEN
1726 5868270 : r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1727 586827 : d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1728 586827 : rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1729 : ELSE
1730 3812310 : r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1731 381231 : d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1732 381231 : rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1733 : END IF
1734 968058 : ka1 = lindh_angle_factor*rho12*rho23; ka2 = lindh_angle_factor*rho23*rho31; ka3 = lindh_angle_factor*rho31*rho12
1735 9680580 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1736 968058 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1737 968058 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1738 968058 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), lindh_denom_floor)
1739 968058 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), lindh_denom_floor)
1740 968058 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), lindh_denom_floor)
1741 968058 : D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1742 968058 : D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1743 968058 : D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1744 968058 : D_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
1745 : D_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
1746 968058 : rsst3*r12(idir)/(d31*d12**3)
1747 : D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1748 968058 : rsst3*r12(jdir)/(d31*d12**3)
1749 968058 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1750 968058 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1751 968058 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1752 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1753 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1754 3987909 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1755 :
1756 : END DO
1757 : END DO
1758 : ELSE
1759 6345666 : DO i = 1, natom
1760 6132546 : IF (i == iat_der .OR. i == jat_der) CYCLE
1761 5706306 : IF (jat_der < iat_der) THEN
1762 2853153 : iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1763 : ELSE
1764 2853153 : iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1765 : END IF
1766 5706306 : IF (jat < i .OR. iat > i) THEN
1767 45387900 : r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1768 4538790 : d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1769 4538790 : rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1770 : ELSE
1771 11675160 : r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1772 1167516 : d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1773 1167516 : rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1774 : END IF
1775 5706306 : ka1 = lindh_angle_factor*rho12*rho23; ka2 = lindh_angle_factor*rho23*rho31; ka3 = lindh_angle_factor*rho31*rho12
1776 57063060 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1777 5706306 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1778 5706306 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1779 5706306 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), lindh_denom_floor)
1780 5706306 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), lindh_denom_floor)
1781 5706306 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), lindh_denom_floor)
1782 5706306 : D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1783 5706306 : D_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
1784 : D_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
1785 5706306 : rsst3*r12(idr)/(d31*d12**3)
1786 5706306 : IF (jat < i .OR. iat > i) THEN
1787 : D_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
1788 4538790 : rsst1*r23(jdr)/(d12*d23**3)
1789 4538790 : D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1790 4538790 : D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1791 : ELSE
1792 1167516 : D_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
1793 : D_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
1794 1167516 : rsst2*r31(jdr)/(d23*d31**3)
1795 1167516 : D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1796 : END IF
1797 5706306 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1798 5706306 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1799 5706306 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1800 :
1801 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1802 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1803 6345666 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1804 : END DO
1805 : END IF
1806 242217 : deriv = 0.25_dp*deriv
1807 :
1808 242217 : END FUNCTION angle_second_deriv
1809 :
1810 : END MODULE bfgs_optimizer
|
