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