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