(git:6d276e9)
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_types, ONLY: cell_type,&
22 pbc
23 USE constraint_fxd, ONLY: fix_atom_control
24 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
25 cp_blacs_env_release,&
26 cp_blacs_env_type
27 USE cp_external_control, ONLY: external_control
28 USE cp_files, ONLY: close_file,&
29 open_file
30 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
31 cp_fm_transpose
32 USE cp_fm_diag, ONLY: choose_eigv_solver
33 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 cp_fm_struct_release,&
35 cp_fm_struct_type
36 USE cp_fm_types, ONLY: cp_fm_create,&
37 cp_fm_get_info,&
38 cp_fm_read_unformatted,&
39 cp_fm_release,&
40 cp_fm_set_all,&
41 cp_fm_to_fm,&
42 cp_fm_type,&
43 cp_fm_write_unformatted
44 USE cp_log_handling, ONLY: cp_get_default_logger,&
45 cp_logger_type,&
46 cp_to_string
47 USE cp_output_handling, ONLY: cp_iterate,&
48 cp_p_file,&
49 cp_print_key_finished_output,&
50 cp_print_key_should_output,&
51 cp_print_key_unit_nr
52 USE cp_subsys_types, ONLY: cp_subsys_get,&
53 cp_subsys_type
54 USE force_env_types, ONLY: force_env_get,&
55 force_env_type
56 USE global_types, ONLY: global_environment_type
57 USE gopt_f_methods, ONLY: cp_eval_at,&
58 gopt_f_ii,&
59 gopt_f_io,&
60 gopt_f_io_finalize,&
61 gopt_f_io_init,&
62 print_geo_opt_header,&
63 print_geo_opt_nc
64 USE gopt_f_types, ONLY: gopt_f_type
65 USE gopt_param_types, ONLY: gopt_param_type
66 USE input_constants, ONLY: default_cell_method_id,&
67 default_ts_method_id
68 USE input_section_types, ONLY: section_vals_get_subs_vals,&
69 section_vals_type,&
70 section_vals_val_get,&
71 section_vals_val_set
72 USE kinds, ONLY: default_path_length,&
73 dp
74 USE machine, ONLY: m_flush,&
75 m_walltime
76 USE message_passing, ONLY: mp_para_env_type
77 USE parallel_gemm_api, ONLY: parallel_gemm
78 USE particle_list_types, ONLY: particle_list_type
79 USE space_groups, ONLY: identify_space_group,&
80 print_spgr,&
81 spgr_apply_rotations_coord,&
82 spgr_apply_rotations_force
83 USE space_groups_types, ONLY: spgr_type
84#include "../base/base_uses.f90"
85
86 IMPLICIT NONE
87 PRIVATE
88
89 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bfgs_optimizer'
90 LOGICAL, PARAMETER :: debug_this_module = .true.
91
92 PUBLIC :: geoopt_bfgs
93
94CONTAINS
95
96! **************************************************************************************************
97!> \brief Main driver for BFGS geometry optimizations
98!> \param force_env ...
99!> \param gopt_param ...
100!> \param globenv ...
101!> \param geo_section ...
102!> \param gopt_env ...
103!> \param x0 ...
104!> \par History
105!> 01.2020 modified to perform Space Group Symmetry [pcazade]
106! **************************************************************************************************
107 RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
108
109 TYPE(force_env_type), POINTER :: force_env
110 TYPE(gopt_param_type), POINTER :: gopt_param
111 TYPE(global_environment_type), POINTER :: globenv
112 TYPE(section_vals_type), POINTER :: geo_section
113 TYPE(gopt_f_type), POINTER :: gopt_env
114 REAL(kind=dp), DIMENSION(:), POINTER :: x0
115
116 CHARACTER(len=*), PARAMETER :: routinen = 'geoopt_bfgs'
117 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
118
119 CHARACTER(LEN=5) :: wildcard
120 CHARACTER(LEN=default_path_length) :: hes_filename
121 INTEGER :: handle, hesunit_read, indf, info, &
122 iter_nr, its, maxiter, ndf, nfree, &
123 output_unit
124 LOGICAL :: conv, hesrest, ionode, shell_present, &
125 should_stop, use_mod_hes, use_rfo
126 REAL(kind=dp) :: ediff, emin, eold, etot, pred, rad, rat, &
127 step, t_diff, t_now, t_old
128 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
129 REAL(kind=dp), DIMENSION(:), POINTER :: g
130 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
131 TYPE(cp_blacs_env_type), POINTER :: blacs_env
132 TYPE(cp_fm_struct_type), POINTER :: fm_struct_hes
133 TYPE(cp_fm_type) :: eigvec_mat, hess_mat, hess_tmp
134 TYPE(cp_logger_type), POINTER :: logger
135 TYPE(cp_subsys_type), POINTER :: subsys
136 TYPE(mp_para_env_type), POINTER :: para_env
137 TYPE(section_vals_type), POINTER :: print_key, root_section
138 TYPE(spgr_type), POINTER :: spgr
139
140 NULLIFY (logger, g, blacs_env, spgr)
141 logger => cp_get_default_logger()
142 para_env => force_env%para_env
143 root_section => force_env%root_section
144 spgr => gopt_env%spgr
145 t_old = m_walltime()
146
147 CALL timeset(routinen, handle)
148 CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
149 print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
150 ionode = para_env%is_source()
151 maxiter = gopt_param%max_iter
152 conv = .false.
153 rat = 0.0_dp
154 wildcard = " BFGS"
155 hes_filename = ""
156
157 ! Stop if not yet implemented
158 SELECT CASE (gopt_env%type_id)
159 CASE (default_ts_method_id)
160 cpabort("BFGS method not yet working with DIMER")
161 END SELECT
162
163 CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
164 CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
165 CALL section_vals_val_get(geo_section, "BFGS%RESTART_HESSIAN", l_val=hesrest)
166 output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
167 extension=".geoLog")
168 IF (output_unit > 0) THEN
169 IF (use_rfo) THEN
170 WRITE (unit=output_unit, fmt="(/,T2,A,T78,A3)") &
171 "BFGS| Use rational function optimization for step estimation: ", "YES"
172 ELSE
173 WRITE (unit=output_unit, fmt="(/,T2,A,T78,A3)") &
174 "BFGS| Use rational function optimization for step estimation: ", " NO"
175 END IF
176 IF (use_mod_hes) THEN
177 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
178 "BFGS| Use model Hessian for initial guess: ", "YES"
179 ELSE
180 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
181 "BFGS| Use model Hessian for initial guess: ", " NO"
182 END IF
183 IF (hesrest) THEN
184 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
185 "BFGS| Restart Hessian: ", "YES"
186 ELSE
187 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
188 "BFGS| Restart Hessian: ", " NO"
189 END IF
190 WRITE (unit=output_unit, fmt="(T2,A,T61,F20.3)") &
191 "BFGS| Trust radius: ", rad
192 END IF
193
194 ndf = SIZE(x0)
195 nfree = gopt_env%nfree
196 IF (ndf > 3000) THEN
197 CALL cp_warn(__location__, &
198 "The dimension of the Hessian matrix ("// &
199 trim(adjustl(cp_to_string(ndf)))//") is greater than 3000. "// &
200 "The diagonalisation of the full Hessian matrix needed for BFGS "// &
201 "is computationally expensive. You should consider to use the linear "// &
202 "scaling variant L-BFGS instead.")
203 END IF
204
205 ! Initialize hessian (hes = unitary matrix or model hessian )
206 CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
207 globenv%blacs_repeatable)
208 CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
209 nrow_global=ndf, ncol_global=ndf)
210 CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
211 CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
212 CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
213 ALLOCATE (eigval(ndf))
214 eigval(:) = 0.0_dp
215
216 CALL force_env_get(force_env=force_env, subsys=subsys)
217 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
218 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
219 IF (use_mod_hes) THEN
220 IF (shell_present) THEN
221 CALL cp_warn(__location__, &
222 "No model Hessian is available for core-shell models. "// &
223 "A unit matrix is used as the initial Hessian.")
224 use_mod_hes = .false.
225 END IF
226 IF (gopt_env%type_id == default_cell_method_id) THEN
227 CALL cp_warn(__location__, &
228 "No model Hessian is available for cell optimizations. "// &
229 "A unit matrix is used as the initial Hessian.")
230 use_mod_hes = .false.
231 END IF
232 END IF
233
234 IF (use_mod_hes) THEN
235 CALL cp_fm_set_all(hess_mat, alpha=zero)
236 CALL construct_initial_hess(gopt_env%force_env, hess_mat)
237 CALL cp_fm_to_fm(hess_mat, hess_tmp)
238 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info)
239 ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?)
240 IF (info /= 0) THEN
241 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
242 IF (output_unit > 0) THEN
243 WRITE (output_unit, *) &
244 "BFGS: Matrix diagonalization failed, using unity as model Hessian."
245 END IF
246 ELSE
247 DO its = 1, SIZE(eigval)
248 IF (eigval(its) < 0.1_dp) eigval(its) = 0.1_dp
249 END DO
250 CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
251 CALL cp_fm_column_scale(eigvec_mat, eigval)
252 CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
253 END IF
254 ELSE
255 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
256 END IF
257
258 ALLOCATE (xold(ndf))
259 xold(:) = x0(:)
260
261 ALLOCATE (g(ndf))
262 g(:) = 0.0_dp
263
264 ALLOCATE (gold(ndf))
265 gold(:) = 0.0_dp
266
267 ALLOCATE (dx(ndf))
268 dx(:) = 0.0_dp
269
270 ALLOCATE (dg(ndf))
271 dg(:) = 0.0_dp
272
273 ALLOCATE (work(ndf))
274 work(:) = 0.0_dp
275
276 ALLOCATE (dr(ndf))
277 dr(:) = 0.0_dp
278
279 ! find space_group
280 CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
281 IF (spgr%keep_space_group) THEN
282 CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
283 CALL spgr_apply_rotations_coord(spgr, x0)
284 CALL print_spgr(spgr)
285 END IF
286
287 ! Geometry optimization starts now
288 CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
289 CALL print_geo_opt_header(gopt_env, output_unit, wildcard)
290
291 ! Calculate Energy & Gradients
292 CALL cp_eval_at(gopt_env, x0, etot, g, master=gopt_env%force_env%para_env%mepos, &
293 para_env=gopt_env%force_env%para_env)
294
295 ! Symmetrize coordinates and forces
296 IF (spgr%keep_space_group) THEN
297 CALL spgr_apply_rotations_coord(spgr, x0)
298 CALL spgr_apply_rotations_force(spgr, g)
299 END IF
300
301 ! Print info at time 0
302 emin = etot
303 t_now = m_walltime()
304 t_diff = t_now - t_old
305 t_old = t_now
306 CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
307 DO its = iter_nr + 1, maxiter
308 CALL cp_iterate(logger%iter_info, last=(its == maxiter))
309 CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
310 CALL gopt_f_ii(its, output_unit)
311
312 ! Hessian update/restarting
313 IF (((its - iter_nr) == 1) .AND. hesrest) THEN
314 IF (ionode) THEN
315 CALL section_vals_val_get(geo_section, "BFGS%RESTART_FILE_NAME", c_val=hes_filename)
316 IF (len_trim(hes_filename) == 0) THEN
317 ! Set default Hessian restart file name if no file name is defined
318 hes_filename = trim(logger%iter_info%project_name)//"-BFGS.Hessian"
319 END IF
320 IF (output_unit > 0) THEN
321 WRITE (unit=output_unit, fmt="(/,T2,A)") &
322 "BFGS| Checking for Hessian restart file <"//trim(adjustl(hes_filename))//">"
323 END IF
324 CALL open_file(file_name=trim(hes_filename), file_status="OLD", &
325 file_form="UNFORMATTED", file_action="READ", unit_number=hesunit_read)
326 IF (output_unit > 0) THEN
327 WRITE (unit=output_unit, fmt="(T2,A)") &
328 "BFGS| Hessian restart file read"
329 END IF
330 END IF
331 CALL cp_fm_read_unformatted(hess_mat, hesunit_read)
332 IF (ionode) CALL close_file(unit_number=hesunit_read)
333 ELSE
334 IF ((its - iter_nr) > 1) THEN
335 ! Symmetrize old coordinates and old forces
336 IF (spgr%keep_space_group) THEN
337 CALL spgr_apply_rotations_coord(spgr, xold)
338 CALL spgr_apply_rotations_force(spgr, gold)
339 END IF
340
341 DO indf = 1, ndf
342 dx(indf) = x0(indf) - xold(indf)
343 dg(indf) = g(indf) - gold(indf)
344 END DO
345
346 CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
347
348 ! Symmetrize coordinates and forces change
349 IF (spgr%keep_space_group) THEN
350 CALL spgr_apply_rotations_force(spgr, dx)
351 CALL spgr_apply_rotations_force(spgr, dg)
352 END IF
353
354 !Possibly dump the Hessian file
355 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
356 CALL write_bfgs_hessian(geo_section, hess_mat, logger)
357 END IF
358 END IF
359 END IF
360
361 ! Symmetrize coordinates and forces
362 IF (spgr%keep_space_group) THEN
363 CALL spgr_apply_rotations_coord(spgr, x0)
364 CALL spgr_apply_rotations_force(spgr, g)
365 END IF
366
367 ! Setting the present positions & gradients as old
368 xold(:) = x0
369 gold(:) = g
370
371 ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
372 CALL cp_fm_to_fm(hess_mat, hess_tmp)
373
374 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval, info=info)
375
376 ! In rare cases the diagonalization of hess_mat fails (bug in scalapack?)
377 IF (info /= 0) THEN
378 IF (output_unit > 0) THEN
379 WRITE (output_unit, *) &
380 "BFGS: Matrix diagonalization failed, resetting Hessian to unity."
381 END IF
382 CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
383 CALL cp_fm_to_fm(hess_mat, hess_tmp)
384 CALL choose_eigv_solver(hess_tmp, eigvec_mat, eigval)
385 END IF
386
387 IF (use_rfo) THEN
388 CALL set_hes_eig(ndf, eigval, work)
389 dx(:) = eigval
390 CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
391 END IF
392 CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
393
394 ! Symmetrize dr
395 IF (spgr%keep_space_group) THEN
396 CALL spgr_apply_rotations_force(spgr, dr)
397 END IF
398
399 CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
400
401 ! Update the atomic positions
402 x0 = x0 + dr
403
404 ! Symmetrize coordinates
405 IF (spgr%keep_space_group) THEN
406 CALL spgr_apply_rotations_coord(spgr, x0)
407 END IF
408
409 CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
410 eold = etot
411
412 ! Energy & Gradients at new step
413 CALL cp_eval_at(gopt_env, x0, etot, g, master=gopt_env%force_env%para_env%mepos, &
414 para_env=gopt_env%force_env%para_env)
415
416 ediff = etot - eold
417
418 ! Symmetrize forces
419 IF (spgr%keep_space_group) THEN
420 CALL spgr_apply_rotations_force(spgr, g)
421 END IF
422
423 ! check for an external exit command
424 CALL external_control(should_stop, "GEO", globenv=globenv)
425 IF (should_stop) EXIT
426
427 ! Some IO and Convergence check
428 t_now = m_walltime()
429 t_diff = t_now - t_old
430 t_old = t_now
431 CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, &
432 eold, emin, wildcard, gopt_param, ndf, dr, g, conv, pred, rat, &
433 step, rad, used_time=t_diff)
434
435 IF (conv .OR. (its == maxiter)) EXIT
436 IF (etot < emin) emin = etot
437 IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
438 END DO
439
440 IF (its == maxiter .AND. (.NOT. conv)) THEN
441 CALL print_geo_opt_nc(gopt_env, output_unit)
442 END IF
443
444 ! show space_group
445 CALL section_vals_val_get(geo_section, "SHOW_SPACE_GROUP", l_val=spgr%show_space_group)
446 IF (spgr%show_space_group) THEN
447 CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
448 CALL print_spgr(spgr)
449 END IF
450
451 ! Write final information, if converged
452 CALL cp_iterate(logger%iter_info, last=.true., increment=0)
453 CALL write_bfgs_hessian(geo_section, hess_mat, logger)
454 CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
455 gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
456
457 CALL cp_fm_struct_release(fm_struct_hes)
458 CALL cp_fm_release(hess_mat)
459 CALL cp_fm_release(eigvec_mat)
460 CALL cp_fm_release(hess_tmp)
461
462 CALL cp_blacs_env_release(blacs_env)
463 DEALLOCATE (xold)
464 DEALLOCATE (g)
465 DEALLOCATE (gold)
466 DEALLOCATE (dx)
467 DEALLOCATE (dg)
468 DEALLOCATE (eigval)
469 DEALLOCATE (work)
470 DEALLOCATE (dr)
471
472 CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
473 "PRINT%PROGRAM_RUN_INFO")
474 CALL timestop(handle)
475
476 END SUBROUTINE geoopt_bfgs
477
478! **************************************************************************************************
479!> \brief ...
480!> \param ndf ...
481!> \param dg ...
482!> \param eigval ...
483!> \param work ...
484!> \param eigvec_mat ...
485!> \param g ...
486!> \param para_env ...
487! **************************************************************************************************
488 SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
489
490 INTEGER, INTENT(IN) :: ndf
491 REAL(kind=dp), INTENT(INOUT) :: dg(ndf), eigval(ndf), work(ndf)
492 TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat
493 REAL(kind=dp), INTENT(INOUT) :: g(ndf)
494 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
495
496 CHARACTER(LEN=*), PARAMETER :: routinen = 'rat_fun_opt'
497 REAL(kind=dp), PARAMETER :: one = 1.0_dp
498
499 INTEGER :: handle, i, indf, iref, iter, j, k, l, &
500 maxit, ncol_local, nrow_local
501 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
502 LOGICAL :: bisec, conv, fail, set
503 REAL(kind=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
504 ln, lp, ssize, step, stol
505 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
506 POINTER :: local_data
507
508 CALL timeset(routinen, handle)
509
510 stol = 1.0e-8_dp
511 ssize = 0.2_dp
512 maxit = 999
513 fail = .false.
514 bisec = .false.
515
516 dg = 0._dp
517
518 CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
519 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
520
521 DO i = 1, nrow_local
522 j = row_indices(i)
523 DO k = 1, ncol_local
524 l = col_indices(k)
525 dg(l) = dg(l) + local_data(i, k)*g(j)
526 END DO
527 END DO
528 CALL para_env%sum(dg)
529
530 set = .false.
531
532 DO
533
534! calculating Lambda
535
536 lp = 0.0_dp
537 iref = 1
538 ln = 0.0_dp
539 IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
540
541 conv = .false.
542 iter = 0
543 DO
544 iter = iter + 1
545 fun = 0.0_dp
546 fung = 0.0_dp
547 DO indf = 1, ndf
548 fun = fun + dg(indf)**2/(ln - eigval(indf))
549 fung = fung - dg(indf)**2/((ln - eigval(indf))**2)
550 END DO
551 fun = fun - ln
552 fung = fung - one
553 step = fun/fung
554 ln = ln - step
555 IF (abs(step) < stol) THEN
556 conv = .true.
557 EXIT
558 END IF
559 IF (iter >= maxit) EXIT
560 END DO
561 outer: DO
562 IF (.NOT. conv) THEN
563 conv = .false.
564 bisec = .true.
565 iter = 0
566 maxit = 9999
567 lam1 = 0.0_dp
568 IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
569 fun1 = 0.0_dp
570 DO indf = 1, ndf
571 fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
572 END DO
573 fun1 = fun1 - lam1
574 step = abs(lam1)/1000.0_dp
575 IF (step < ssize) step = ssize
576 inner: DO
577 iter = iter + 1
578 IF (iter > maxit) THEN
579 ln = 0.0_dp
580 lp = 0.0_dp
581 fail = .true.
582 EXIT outer
583 END IF
584 fun2 = 0.0_dp
585 lam2 = lam1 - iter*step
586 DO indf = 1, ndf
587 fun2 = fun2 + dg(indf)**2/(lam2 - eigval(indf))
588 END DO
589 fun2 = fun2 - lam2
590 IF (fun2*fun1 < 0.0_dp) THEN
591 iter = 0
592 DO
593 iter = iter + 1
594 IF (iter > maxit) THEN
595 ln = 0.0_dp
596 lp = 0.0_dp
597 fail = .true.
598 EXIT outer
599 END IF
600 step = (lam1 + lam2)/2
601 fun3 = 0.0_dp
602 DO indf = 1, ndf
603 fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
604 END DO
605 fun3 = fun3 - step
606
607 IF (abs(step - lam2) < stol) THEN
608 ln = step
609 EXIT inner
610 END IF
611
612 IF (fun3*fun1 < stol) THEN
613 lam2 = step
614 ELSE
615 lam1 = step
616 END IF
617 END DO
618 END IF
619 END DO inner
620 END IF
621 IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
622 (eigval(iref) > 0.0_dp))) THEN
623
624 IF (.NOT. bisec) THEN
625 conv = .false.
626 cycle outer
627 END IF
628 ln = 0.0_dp
629 lp = 0.0_dp
630 fail = .true.
631 END IF
632 EXIT outer
633 END DO outer
634
635 IF (fail .AND. .NOT. set) THEN
636 set = .true.
637 DO indf = 1, ndf
638 eigval(indf) = eigval(indf)*work(indf)
639 END DO
640 cycle
641 END IF
642
643 IF (.NOT. set) THEN
644 work(1:ndf) = one
645 END IF
646
647 DO indf = 1, ndf
648 eigval(indf) = eigval(indf) - ln
649 END DO
650 EXIT
651 END DO
652
653 CALL timestop(handle)
654
655 END SUBROUTINE rat_fun_opt
656
657! **************************************************************************************************
658!> \brief ...
659!> \param ndf ...
660!> \param dx ...
661!> \param dg ...
662!> \param hess_mat ...
663!> \param work ...
664!> \param para_env ...
665! **************************************************************************************************
666 SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env)
667 INTEGER, INTENT(IN) :: ndf
668 REAL(kind=dp), INTENT(INOUT) :: dx(ndf), dg(ndf)
669 TYPE(cp_fm_type), INTENT(IN) :: hess_mat
670 REAL(kind=dp), INTENT(INOUT) :: work(ndf)
671 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
672
673 CHARACTER(LEN=*), PARAMETER :: routinen = 'bfgs'
674 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
675
676 INTEGER :: handle, i, j, k, l, ncol_local, &
677 nrow_local
678 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
679 REAL(kind=dp) :: ddot, dxw, gdx
680 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
681 POINTER :: local_hes
682
683 CALL timeset(routinen, handle)
684
685 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
686 local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
687
688 work = zero
689 DO i = 1, nrow_local
690 j = row_indices(i)
691 DO k = 1, ncol_local
692 l = col_indices(k)
693 work(j) = work(j) + local_hes(i, k)*dx(l)
694 END DO
695 END DO
696
697 CALL para_env%sum(work)
698
699 gdx = ddot(ndf, dg, 1, dx, 1)
700 gdx = one/gdx
701 dxw = ddot(ndf, dx, 1, work, 1)
702 dxw = one/dxw
703
704 DO i = 1, nrow_local
705 j = row_indices(i)
706 DO k = 1, ncol_local
707 l = col_indices(k)
708 local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
709 dxw*work(j)*work(l)
710 END DO
711 END DO
712
713 CALL timestop(handle)
714
715 END SUBROUTINE bfgs
716
717! **************************************************************************************************
718!> \brief ...
719!> \param ndf ...
720!> \param eigval ...
721!> \param work ...
722! **************************************************************************************************
723 SUBROUTINE set_hes_eig(ndf, eigval, work)
724 INTEGER, INTENT(IN) :: ndf
725 REAL(kind=dp), INTENT(INOUT) :: eigval(ndf), work(ndf)
726
727 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_hes_eig'
728 REAL(kind=dp), PARAMETER :: max_neg = -0.5_dp, max_pos = 5.0_dp, &
729 min_eig = 0.005_dp, one = 1.0_dp
730
731 INTEGER :: handle, indf
732 LOGICAL :: neg
733
734 CALL timeset(routinen, handle)
735
736 DO indf = 1, ndf
737 IF (eigval(indf) < 0.0_dp) neg = .true.
738 IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
739 END DO
740 DO indf = 1, ndf
741 IF (eigval(indf) < 0.0_dp) THEN
742 IF (eigval(indf) < max_neg) THEN
743 eigval(indf) = max_neg
744 ELSE IF (eigval(indf) > -min_eig) THEN
745 eigval(indf) = -min_eig
746 END IF
747 ELSE IF (eigval(indf) < 1000.0_dp) THEN
748 IF (eigval(indf) < min_eig) THEN
749 eigval(indf) = min_eig
750 ELSE IF (eigval(indf) > max_pos) THEN
751 eigval(indf) = max_pos
752 END IF
753 END IF
754 END DO
755
756 DO indf = 1, ndf
757 IF (eigval(indf) < 0.0_dp) THEN
758 work(indf) = -one
759 ELSE
760 work(indf) = one
761 END IF
762 END DO
763
764 CALL timestop(handle)
765
766 END SUBROUTINE set_hes_eig
767
768! **************************************************************************************************
769!> \brief ...
770!> \param ndf ...
771!> \param eigval ...
772!> \param eigvec_mat ...
773!> \param hess_tmp ...
774!> \param dr ...
775!> \param g ...
776!> \param para_env ...
777!> \param use_rfo ...
778! **************************************************************************************************
779 SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
780
781 INTEGER, INTENT(IN) :: ndf
782 REAL(kind=dp), INTENT(INOUT) :: eigval(ndf)
783 TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat, hess_tmp
784 REAL(kind=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
785 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
786 LOGICAL :: use_rfo
787
788 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
789
790 INTEGER :: i, indf, j, k, l, ncol_local, nrow_local
791 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
792 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
793 POINTER :: local_data
794 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
795 TYPE(cp_fm_type) :: tmp
796
797 CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
798 IF (use_rfo) THEN
799 DO indf = 1, ndf
800 eigval(indf) = one/eigval(indf)
801 END DO
802 ELSE
803 DO indf = 1, ndf
804 eigval(indf) = one/max(0.0001_dp, eigval(indf))
805 END DO
806 END IF
807
808 CALL cp_fm_column_scale(hess_tmp, eigval)
809 CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
810 CALL cp_fm_create(tmp, matrix_struct, name="tmp")
811 CALL cp_fm_set_all(tmp, alpha=zero)
812
813 CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
814
815 CALL cp_fm_transpose(tmp, hess_tmp)
816 CALL cp_fm_release(tmp)
817
818 ! New step
819
820 CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
821 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
822
823 dr = 0.0_dp
824 DO i = 1, nrow_local
825 j = row_indices(i)
826 DO k = 1, ncol_local
827 l = col_indices(k)
828 dr(j) = dr(j) - local_data(i, k)*g(l)
829 END DO
830 END DO
831
832 CALL para_env%sum(dr)
833
834 END SUBROUTINE geoopt_get_step
835
836! **************************************************************************************************
837!> \brief ...
838!> \param ndf ...
839!> \param step ...
840!> \param rad ...
841!> \param rat ...
842!> \param dr ...
843!> \param output_unit ...
844! **************************************************************************************************
845 SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
846 INTEGER, INTENT(IN) :: ndf
847 REAL(kind=dp), INTENT(INOUT) :: step, rad, rat, dr(ndf)
848 INTEGER, INTENT(IN) :: output_unit
849
850 CHARACTER(LEN=*), PARAMETER :: routinen = 'trust_radius'
851 REAL(kind=dp), PARAMETER :: one = 1.0_dp
852
853 INTEGER :: handle
854 REAL(kind=dp) :: scal
855
856 CALL timeset(routinen, handle)
857
858 step = maxval(abs(dr))
859 scal = max(one, rad/step)
860
861 IF (step > rad) THEN
862 rat = rad/step
863 CALL dscal(ndf, rat, dr, 1)
864 step = rad
865 IF (output_unit > 0) THEN
866 WRITE (unit=output_unit, fmt="(/,T2,A,F8.5)") &
867 " Step is scaled; Scaling factor = ", rat
868 CALL m_flush(output_unit)
869 END IF
870 END IF
871 CALL timestop(handle)
872
873 END SUBROUTINE trust_radius
874
875! **************************************************************************************************
876!> \brief ...
877!> \param ndf ...
878!> \param work ...
879!> \param hess_mat ...
880!> \param dr ...
881!> \param g ...
882!> \param conv ...
883!> \param pred ...
884!> \param para_env ...
885! **************************************************************************************************
886 SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
887
888 INTEGER, INTENT(IN) :: ndf
889 REAL(kind=dp), INTENT(INOUT) :: work(ndf)
890 TYPE(cp_fm_type), INTENT(IN) :: hess_mat
891 REAL(kind=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
892 LOGICAL, INTENT(INOUT) :: conv
893 REAL(kind=dp), INTENT(INOUT) :: pred
894 TYPE(mp_para_env_type), POINTER :: para_env
895
896 CHARACTER(LEN=*), PARAMETER :: routinen = 'energy_predict'
897 REAL(kind=dp), PARAMETER :: zero = 0.0_dp
898
899 INTEGER :: handle, i, j, k, l, ncol_local, &
900 nrow_local
901 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
902 REAL(kind=dp) :: ddot, ener1, ener2
903 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
904 POINTER :: local_data
905
906 CALL timeset(routinen, handle)
907
908 ener1 = ddot(ndf, g, 1, dr, 1)
909
910 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
911 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
912
913 work = zero
914 DO i = 1, nrow_local
915 j = row_indices(i)
916 DO k = 1, ncol_local
917 l = col_indices(k)
918 work(j) = work(j) + local_data(i, k)*dr(l)
919 END DO
920 END DO
921
922 CALL para_env%sum(work)
923 ener2 = ddot(ndf, dr, 1, work, 1)
924 pred = ener1 + 0.5_dp*ener2
925 conv = .false.
926 CALL timestop(handle)
927
928 END SUBROUTINE energy_predict
929
930! **************************************************************************************************
931!> \brief ...
932!> \param rat ...
933!> \param rad ...
934!> \param step ...
935!> \param ediff ...
936! **************************************************************************************************
937 SUBROUTINE update_trust_rad(rat, rad, step, ediff)
938
939 REAL(kind=dp), INTENT(INOUT) :: rat, rad, step, ediff
940
941 CHARACTER(LEN=*), PARAMETER :: routinen = 'update_trust_rad'
942 REAL(kind=dp), PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp
943
944 INTEGER :: handle
945
946 CALL timeset(routinen, handle)
947
948 IF (rat > 4.0_dp) THEN
949 IF (ediff < 0.0_dp) THEN
950 rad = step*0.5_dp
951 ELSE
952 rad = step*0.25_dp
953 END IF
954 ELSE IF (rat > 2.0_dp) THEN
955 IF (ediff < 0.0_dp) THEN
956 rad = step*0.75_dp
957 ELSE
958 rad = step*0.5_dp
959 END IF
960 ELSE IF (rat > 4.0_dp/3.0_dp) THEN
961 IF (ediff < 0.0_dp) THEN
962 rad = step
963 ELSE
964 rad = step*0.75_dp
965 END IF
966 ELSE IF (rat > 10.0_dp/9.0_dp) THEN
967 IF (ediff < 0.0_dp) THEN
968 rad = step*1.25_dp
969 ELSE
970 rad = step
971 END IF
972 ELSE IF (rat > 0.9_dp) THEN
973 IF (ediff < 0.0_dp) THEN
974 rad = step*1.5_dp
975 ELSE
976 rad = step*1.25_dp
977 END IF
978 ELSE IF (rat > 0.75_dp) THEN
979 IF (ediff < 0.0_dp) THEN
980 rad = step*1.25_dp
981 ELSE
982 rad = step
983 END IF
984 ELSE IF (rat > 0.5_dp) THEN
985 IF (ediff < 0.0_dp) THEN
986 rad = step
987 ELSE
988 rad = step*0.75_dp
989 END IF
990 ELSE IF (rat > 0.25_dp) THEN
991 IF (ediff < 0.0_dp) THEN
992 rad = step*0.75_dp
993 ELSE
994 rad = step*0.5_dp
995 END IF
996 ELSE IF (ediff < 0.0_dp) THEN
997 rad = step*0.5_dp
998 ELSE
999 rad = step*0.25_dp
1000 END IF
1001
1002 rad = max(rad, min_trust)
1003 rad = min(rad, max_trust)
1004 CALL timestop(handle)
1005
1006 END SUBROUTINE update_trust_rad
1007
1008! **************************************************************************************************
1009
1010! **************************************************************************************************
1011!> \brief ...
1012!> \param geo_section ...
1013!> \param hess_mat ...
1014!> \param logger ...
1015! **************************************************************************************************
1016 SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger)
1017
1018 TYPE(section_vals_type), POINTER :: geo_section
1019 TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1020 TYPE(cp_logger_type), POINTER :: logger
1021
1022 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_bfgs_hessian'
1023
1024 INTEGER :: handle, hesunit
1025
1026 CALL timeset(routinen, handle)
1027
1028 hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", &
1029 extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", &
1030 file_position="REWIND")
1031
1032 CALL cp_fm_write_unformatted(hess_mat, hesunit)
1033
1034 CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
1035
1036 CALL timestop(handle)
1037
1038 END SUBROUTINE write_bfgs_hessian
1039
1040! **************************************************************************************************
1041!> \brief Constructs model Hessian as described in https://doi.org/10.1016/0009-2614(95)00646-L.
1042!> \param force_env ...
1043!> \param hess_mat ...
1044!> \author Florian Schiffmann
1045! **************************************************************************************************
1046 SUBROUTINE construct_initial_hess(force_env, hess_mat)
1047
1048 TYPE(force_env_type), POINTER :: force_env
1049 TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1050
1051 INTEGER :: i, iat_col, iat_row, iglobal, iind, j, &
1052 jat_row, jglobal, jind, k, natom, &
1053 ncol_local, nrow_local, z
1054 INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row
1055 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1056 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij
1057 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij
1058 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1059 POINTER :: fixed, local_data
1060 REAL(kind=dp), DIMENSION(3, 3) :: alpha, r0
1061 TYPE(cell_type), POINTER :: cell
1062 TYPE(cp_subsys_type), POINTER :: subsys
1063 TYPE(particle_list_type), POINTER :: particles
1064
1065 CALL cite_reference(lindh1995)
1066
1067 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1068 CALL cp_subsys_get(subsys, &
1069 particles=particles)
1070
1071 alpha(1, :) = [1._dp, 0.3949_dp, 0.3949_dp]
1072 alpha(2, :) = [0.3494_dp, 0.2800_dp, 0.2800_dp]
1073 alpha(3, :) = [0.3494_dp, 0.2800_dp, 0.1800_dp]
1074
1075 r0(1, :) = [1.35_dp, 2.10_dp, 2.53_dp]
1076 r0(2, :) = [2.10_dp, 2.87_dp, 3.40_dp]
1077 r0(3, :) = [2.53_dp, 3.40_dp, 3.40_dp]
1078
1079 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1080 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1081 natom = particles%n_els
1082 ALLOCATE (at_row(natom))
1083 ALLOCATE (rho_ij(natom, natom))
1084 ALLOCATE (d_ij(natom, natom))
1085 ALLOCATE (r_ij(natom, natom, 3))
1086 ALLOCATE (fixed(3, natom))
1087 fixed = 1.0_dp
1088 CALL fix_atom_control(force_env, fixed)
1089 DO i = 1, 3
1090 CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1091 END DO
1092 rho_ij = 0
1093 !XXXX insert proper rows !XXX
1094 at_row = 3
1095 DO i = 1, natom
1096 CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
1097 IF (z <= 10) at_row(i) = 2
1098 IF (z <= 2) at_row(i) = 1
1099 END DO
1100 DO i = 2, natom
1101 iat_row = at_row(i)
1102 DO j = 1, i - 1
1103 jat_row = at_row(j)
1104 !pbc for a distance vector
1105 r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
1106 r_ij(i, j, :) = -r_ij(j, i, :)
1107 d_ij(j, i) = norm2(r_ij(j, i, :))
1108 d_ij(i, j) = d_ij(j, i)
1109 rho_ij(j, i) = exp(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1110 rho_ij(i, j) = rho_ij(j, i)
1111 END DO
1112 END DO
1113 DO i = 1, ncol_local
1114 iglobal = col_indices(i)
1115 iind = mod(iglobal - 1, 3) + 1
1116 iat_col = (iglobal + 2)/3
1117 IF (iat_col > natom) cycle
1118 DO j = 1, nrow_local
1119 jglobal = row_indices(j)
1120 jind = mod(jglobal - 1, 3) + 1
1121 iat_row = (jglobal + 2)/3
1122 IF (iat_row > natom) cycle
1123 IF (iat_row /= iat_col) THEN
1124 IF (d_ij(iat_row, iat_col) < 6.0_dp) THEN
1125 local_data(j, i) = local_data(j, i) + &
1126 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1127 END IF
1128 ELSE
1129 local_data(j, i) = local_data(j, i) + &
1130 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1131 END IF
1132 IF (iat_col /= iat_row) THEN
1133 IF (d_ij(iat_row, iat_col) < 6.0_dp) THEN
1134 local_data(j, i) = local_data(j, i) - &
1135 dist_second_deriv(r_ij(iat_col, iat_row, :), &
1136 iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1137 END IF
1138 ELSE
1139 DO k = 1, natom
1140 IF (k == iat_col) cycle
1141 IF (d_ij(iat_row, k) < 6.0_dp) THEN
1142 local_data(j, i) = local_data(j, i) + &
1143 dist_second_deriv(r_ij(iat_col, k, :), &
1144 iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1145 END IF
1146 END DO
1147 END IF
1148 IF (fixed(jind, iat_row) < 0.5_dp .OR. fixed(iind, iat_col) < 0.5_dp) THEN
1149 local_data(j, i) = 0.0_dp
1150 IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
1151 END IF
1152 END DO
1153 END DO
1154 DEALLOCATE (fixed)
1155 DEALLOCATE (rho_ij)
1156 DEALLOCATE (d_ij)
1157 DEALLOCATE (r_ij)
1158 DEALLOCATE (at_row)
1159
1160 END SUBROUTINE construct_initial_hess
1161
1162! **************************************************************************************************
1163!> \brief ...
1164!> \param r1 ...
1165!> \param i ...
1166!> \param j ...
1167!> \param d ...
1168!> \param rho ...
1169!> \return ...
1170! **************************************************************************************************
1171 PURE FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv)
1172 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r1
1173 INTEGER, INTENT(IN) :: i, j
1174 REAL(kind=dp), INTENT(IN) :: d, rho
1175 REAL(kind=dp) :: deriv
1176
1177 deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
1178 END FUNCTION dist_second_deriv
1179
1180! **************************************************************************************************
1181!> \brief ...
1182!> \param r_ij ...
1183!> \param d_ij ...
1184!> \param rho_ij ...
1185!> \param idir ...
1186!> \param jdir ...
1187!> \param iat_der ...
1188!> \param jat_der ...
1189!> \param natom ...
1190!> \return ...
1191! **************************************************************************************************
1192 FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv)
1193 REAL(kind=dp), DIMENSION(:, :, :) :: r_ij
1194 REAL(kind=dp), DIMENSION(:, :) :: d_ij, rho_ij
1195 INTEGER :: idir, jdir, iat_der, jat_der, natom
1196 REAL(kind=dp) :: deriv
1197
1198 INTEGER :: i, iat, idr, j, jat, jdr
1199 REAL(kind=dp) :: d12, d23, d31, d_mat(3, 2), denom1, &
1200 denom2, denom3, ka1, ka2, ka3, rho12, &
1201 rho23, rho31, rsst1, rsst2, rsst3
1202 REAL(kind=dp), DIMENSION(3) :: r12, r23, r31
1203
1204 deriv = 0._dp
1205 IF (iat_der == jat_der) THEN
1206 DO i = 1, natom - 1
1207 IF (rho_ij(iat_der, i) < 0.00001) cycle
1208 DO j = i + 1, natom
1209 IF (rho_ij(iat_der, j) < 0.00001) cycle
1210 IF (i == iat_der .OR. j == iat_der) cycle
1211 IF (iat_der < i .OR. iat_der > j) THEN
1212 r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1213 d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1214 rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1215 ELSE
1216 r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1217 d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1218 rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1219 END IF
1220 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1221 rsst1 = dot_product(r12, r23); rsst2 = dot_product(r23, r31); rsst3 = dot_product(r31, r12)
1222 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1223 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1224 denom1 = sign(1.0_dp, denom1)*max(abs(denom1), 0.01_dp)
1225 denom2 = sign(1.0_dp, denom2)*max(abs(denom2), 0.01_dp)
1226 denom3 = sign(1.0_dp, denom3)*max(abs(denom3), 0.01_dp)
1227 d_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1228 d_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1229 d_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1230 d_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
1231 d_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
1232 rsst3*r12(idir)/(d31*d12**3)
1233 d_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1234 rsst3*r12(jdir)/(d31*d12**3)
1235 IF (abs(denom1) <= 0.011_dp) d_mat(1, 1) = 0.0_dp
1236 IF (abs(denom2) <= 0.011_dp) d_mat(2, 1) = 0.0_dp
1237 IF (abs(denom3) <= 0.011_dp) d_mat(3, 1) = 0.0_dp
1238 deriv = deriv + ka1*d_mat(1, 1)*d_mat(1, 2)/denom1 + &
1239 ka2*d_mat(2, 1)*d_mat(2, 2)/denom2 + &
1240 ka3*d_mat(3, 1)*d_mat(3, 2)/denom3
1241
1242 END DO
1243 END DO
1244 ELSE
1245 DO i = 1, natom
1246 IF (i == iat_der .OR. i == jat_der) cycle
1247 IF (jat_der < iat_der) THEN
1248 iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1249 ELSE
1250 iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1251 END IF
1252 IF (jat < i .OR. iat > i) THEN
1253 r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1254 d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1255 rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1256 ELSE
1257 r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1258 d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1259 rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1260 END IF
1261 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1262 rsst1 = dot_product(r12, r23); rsst2 = dot_product(r23, r31); rsst3 = dot_product(r31, r12)
1263 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1264 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1265 denom1 = sign(1.0_dp, denom1)*max(abs(denom1), 0.01_dp)
1266 denom2 = sign(1.0_dp, denom2)*max(abs(denom2), 0.01_dp)
1267 denom3 = sign(1.0_dp, denom3)*max(abs(denom3), 0.01_dp)
1268 d_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1269 d_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
1270 d_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
1271 rsst3*r12(idr)/(d31*d12**3)
1272 IF (jat < i .OR. iat > i) THEN
1273 d_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
1274 rsst1*r23(jdr)/(d12*d23**3)
1275 d_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1276 d_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1277 ELSE
1278 d_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
1279 d_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
1280 rsst2*r31(jdr)/(d23*d31**3)
1281 d_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1282 END IF
1283 IF (abs(denom1) <= 0.011_dp) d_mat(1, 1) = 0.0_dp
1284 IF (abs(denom2) <= 0.011_dp) d_mat(2, 1) = 0.0_dp
1285 IF (abs(denom3) <= 0.011_dp) d_mat(3, 1) = 0.0_dp
1286
1287 deriv = deriv + ka1*d_mat(1, 1)*d_mat(1, 2)/denom1 + &
1288 ka2*d_mat(2, 1)*d_mat(2, 2)/denom2 + &
1289 ka3*d_mat(3, 1)*d_mat(3, 2)/denom3
1290 END DO
1291 END IF
1292 deriv = 0.25_dp*deriv
1293
1294 END FUNCTION angle_second_deriv
1295
1296END MODULE bfgs_optimizer
cell_types::pbc
Definition cell_types.F:104
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
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::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:108
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_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. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
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:1690
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:1532
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:1038
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2132
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2363
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:1211
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:404
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:659
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:505
gopt_f_methods::print_geo_opt_nc
subroutine, public print_geo_opt_nc(gopt_env, output_unit)
...
Definition gopt_f_methods.F:1250
gopt_f_methods::cp_eval_at
subroutine, public cp_eval_at(gopt_env, x, f, gradient, master, final_evaluation, para_env)
evaluete the potential energy and its gradients using an array with same dimension as the particle_se...
Definition gopt_f_methods.F:173
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:381
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:559
input_constants::default_ts_method_id
integer, parameter, public default_ts_method_id
Definition input_constants.F:559
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:124
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:141
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:59
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:767
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