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