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