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