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 876 : 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 876 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
130 876 : 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 876 : NULLIFY (logger, g, blacs_env, spgr)
142 1752 : logger => cp_get_default_logger()
143 876 : para_env => force_env%para_env
144 876 : root_section => force_env%root_section
145 876 : spgr => gopt_env%spgr
146 876 : t_old = m_walltime()
147 :
148 876 : CALL timeset(routineN, handle)
149 876 : CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
150 876 : print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
151 876 : ionode = para_env%is_source()
152 876 : maxiter = gopt_param%max_iter
153 876 : conv = .FALSE.
154 876 : rat = 0.0_dp
155 876 : wildcard = " BFGS"
156 876 : hes_filename = ""
157 :
158 : ! Stop if not yet implemented
159 876 : SELECT CASE (gopt_env%type_id)
160 : CASE (default_ts_method_id)
161 876 : CPABORT("BFGS method not yet working with DIMER")
162 : END SELECT
163 :
164 876 : CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
165 876 : CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
166 876 : 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 876 : extension=".geoLog")
169 876 : IF (output_unit > 0) THEN
170 454 : 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 449 : "BFGS| Use rational function optimization for step estimation: ", " NO"
176 : END IF
177 454 : 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 59 : "BFGS| Use model Hessian for initial guess: ", " NO"
183 : END IF
184 454 : 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 453 : "BFGS| Restart Hessian: ", " NO"
190 : END IF
191 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
192 454 : "BFGS| Trust radius: ", rad
193 : END IF
194 :
195 876 : ndf = SIZE(x0)
196 876 : nfree = gopt_env%nfree
197 876 : 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 876 : globenv%blacs_repeatable)
209 : CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
210 876 : nrow_global=ndf, ncol_global=ndf)
211 876 : CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
212 876 : CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
213 876 : CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
214 2628 : ALLOCATE (eigval(ndf))
215 876 : eigval(:) = 0.0_dp
216 :
217 876 : CALL force_env_get(force_env=force_env, subsys=subsys)
218 876 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
219 876 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
220 876 : IF (use_mod_hes) THEN
221 758 : 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 758 : 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 876 : IF (use_mod_hes) THEN
236 754 : CALL cp_fm_set_all(hess_mat, alpha=zero)
237 754 : CALL construct_initial_hess(gopt_env%force_env, hess_mat)
238 754 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
239 754 : 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 754 : 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 16420 : DO its = 1, SIZE(eigval)
249 16420 : IF (eigval(its) < 0.1_dp) eigval(its) = 0.1_dp
250 : END DO
251 754 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
252 754 : CALL cp_fm_column_scale(eigvec_mat, eigval)
253 754 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
254 : END IF
255 : ELSE
256 122 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
257 : END IF
258 :
259 2628 : ALLOCATE (xold(ndf))
260 23964 : xold(:) = x0(:)
261 :
262 1752 : ALLOCATE (g(ndf))
263 23964 : g(:) = 0.0_dp
264 :
265 1752 : ALLOCATE (gold(ndf))
266 876 : gold(:) = 0.0_dp
267 :
268 2628 : ALLOCATE (dx(ndf))
269 876 : dx(:) = 0.0_dp
270 :
271 2628 : ALLOCATE (dg(ndf))
272 876 : dg(:) = 0.0_dp
273 :
274 2628 : ALLOCATE (work(ndf))
275 876 : work(:) = 0.0_dp
276 :
277 2628 : ALLOCATE (dr(ndf))
278 876 : dr(:) = 0.0_dp
279 :
280 : ! find space_group
281 876 : CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
282 876 : IF (spgr%keep_space_group) THEN
283 8 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
284 8 : CALL spgr_apply_rotations_coord(spgr, x0)
285 8 : CALL print_spgr(spgr)
286 : END IF
287 :
288 : ! Geometry optimization starts now
289 876 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
290 876 : 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 876 : .FALSE., gopt_env%force_env%para_env)
295 :
296 : ! Symmetrize coordinates and forces
297 876 : IF (spgr%keep_space_group) THEN
298 8 : CALL spgr_apply_rotations_coord(spgr, x0)
299 8 : CALL spgr_apply_rotations_force(spgr, g)
300 : END IF
301 :
302 : ! Print info at time 0
303 876 : emin = etot
304 876 : t_now = m_walltime()
305 876 : t_diff = t_now - t_old
306 876 : t_old = t_now
307 876 : CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
308 4138 : DO its = iter_nr + 1, maxiter
309 4128 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
310 4128 : CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
311 4128 : CALL gopt_f_ii(its, output_unit)
312 :
313 : ! Hessian update/restarting
314 4128 : 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 4126 : IF ((its - iter_nr) > 1) THEN
336 : ! Symmetrize old coordinates and old forces
337 3262 : 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 254887 : DO indf = 1, ndf
343 251625 : dx(indf) = x0(indf) - xold(indf)
344 254887 : dg(indf) = g(indf) - gold(indf)
345 : END DO
346 :
347 3262 : CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
348 :
349 : ! Symmetrize coordinates and forces change
350 3262 : 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 3262 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
357 2663 : 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 4128 : IF (spgr%keep_space_group) THEN
364 8 : CALL spgr_apply_rotations_coord(spgr, x0)
365 8 : CALL spgr_apply_rotations_force(spgr, g)
366 : END IF
367 :
368 : ! Setting the present positions & gradients as old
369 278697 : xold(:) = x0
370 278697 : gold(:) = g
371 :
372 : ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
373 4128 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
374 :
375 4128 : 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 4128 : 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 4128 : IF (use_rfo) THEN
389 801 : CALL set_hes_eig(ndf, eigval, work)
390 80871 : dx(:) = eigval
391 801 : CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
392 : END IF
393 4128 : CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
394 :
395 : ! Symmetrize dr
396 4128 : IF (spgr%keep_space_group) THEN
397 8 : CALL spgr_apply_rotations_force(spgr, dr)
398 : END IF
399 :
400 4128 : CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
401 :
402 : ! Update the atomic positions
403 278697 : x0 = x0 + dr
404 :
405 : ! Symmetrize coordinates
406 4128 : IF (spgr%keep_space_group) THEN
407 8 : CALL spgr_apply_rotations_coord(spgr, x0)
408 : END IF
409 :
410 4128 : CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
411 4128 : 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 4128 : .FALSE., gopt_env%force_env%para_env)
416 :
417 4128 : ediff = etot - eold
418 :
419 : ! Symmetrize forces
420 4128 : IF (spgr%keep_space_group) THEN
421 8 : CALL spgr_apply_rotations_force(spgr, g)
422 : END IF
423 :
424 : ! check for an external exit command
425 4128 : CALL external_control(should_stop, "GEO", globenv=globenv)
426 4128 : IF (should_stop) EXIT
427 :
428 : ! Some IO and Convergence check
429 4128 : t_now = m_walltime()
430 4128 : t_diff = t_now - t_old
431 4128 : 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 4128 : step, rad, used_time=t_diff)
435 :
436 4128 : IF (conv .OR. (its == maxiter)) EXIT
437 3262 : IF (etot < emin) emin = etot
438 12394 : IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
439 : END DO
440 :
441 876 : IF (its == maxiter .AND. (.NOT. conv)) THEN
442 605 : CALL print_geo_opt_nc(gopt_env, output_unit)
443 : END IF
444 :
445 : ! show space_group
446 876 : CALL section_vals_val_get(geo_section, "SHOW_SPACE_GROUP", l_val=spgr%show_space_group)
447 876 : 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 876 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
454 876 : 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 876 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
457 :
458 876 : CALL cp_fm_struct_release(fm_struct_hes)
459 876 : CALL cp_fm_release(hess_mat)
460 876 : CALL cp_fm_release(eigvec_mat)
461 876 : CALL cp_fm_release(hess_tmp)
462 :
463 876 : CALL cp_blacs_env_release(blacs_env)
464 876 : DEALLOCATE (xold)
465 876 : DEALLOCATE (g)
466 876 : DEALLOCATE (gold)
467 876 : DEALLOCATE (dx)
468 876 : DEALLOCATE (dg)
469 876 : DEALLOCATE (eigval)
470 876 : DEALLOCATE (work)
471 876 : DEALLOCATE (dr)
472 :
473 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
474 876 : "PRINT%PROGRAM_RUN_INFO")
475 876 : CALL timestop(handle)
476 :
477 5256 : 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 1602 : 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 801 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
503 : LOGICAL :: bisec, fail, set
504 : REAL(KIND=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
505 : ln, lp, ssize, step, stol
506 801 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
507 :
508 801 : CALL timeset(routineN, handle)
509 :
510 801 : stol = 1.0E-8_dp
511 801 : ssize = 0.2_dp
512 801 : maxit = 999
513 801 : fail = .FALSE.
514 801 : bisec = .FALSE.
515 :
516 80871 : dg = 0._dp
517 :
518 : CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
519 801 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
520 :
521 50241 : DO i = 1, nrow_local
522 49440 : j = row_indices(i)
523 11471241 : DO k = 1, ncol_local
524 11421000 : l = col_indices(k)
525 11470440 : dg(l) = dg(l) + local_data(i, k)*g(j)
526 : END DO
527 : END DO
528 801 : CALL para_env%sum(dg)
529 :
530 801 : set = .FALSE.
531 :
532 : DO
533 :
534 : ! calculating Lambda
535 :
536 801 : lp = 0.0_dp
537 801 : iref = 1
538 801 : ln = 0.0_dp
539 801 : IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
540 :
541 801 : iter = 0
542 : DO
543 2502 : iter = iter + 1
544 2502 : fun = 0.0_dp
545 2502 : fung = 0.0_dp
546 239682 : DO indf = 1, ndf
547 237180 : fun = fun + dg(indf)**2/(ln - eigval(indf))
548 239682 : fung = fung - dg(indf)**2/(ln - eigval(indf)**2)
549 : END DO
550 2502 : fun = fun - ln
551 2502 : fung = fung - one
552 2502 : step = fun/fung
553 2502 : ln = ln - step
554 2502 : IF (ABS(step) < stol) GOTO 200
555 1702 : IF (iter >= maxit) EXIT
556 : END DO
557 : 100 CONTINUE
558 31 : bisec = .TRUE.
559 31 : iter = 0
560 31 : maxit = 9999
561 31 : lam1 = 0.0_dp
562 31 : IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
563 : fun1 = 0.0_dp
564 31 : DO indf = 1, ndf
565 31 : fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
566 : END DO
567 1 : fun1 = fun1 - lam1
568 1 : step = ABS(lam1)/1000.0_dp
569 : IF (step < ssize) step = ssize
570 : DO
571 1 : iter = iter + 1
572 1 : IF (iter > maxit) THEN
573 : ln = 0.0_dp
574 801 : lp = 0.0_dp
575 : fail = .TRUE.
576 : GOTO 300
577 : END IF
578 1 : fun2 = 0.0_dp
579 1 : lam2 = lam1 - iter*step
580 31 : DO indf = 1, ndf
581 31 : fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf))
582 : END DO
583 1 : fun2 = fun2 - lam2
584 802 : IF (fun2*fun1 < 0.0_dp) THEN
585 : iter = 0
586 : DO
587 25 : iter = iter + 1
588 25 : IF (iter > maxit) THEN
589 : ln = 0.0_dp
590 : lp = 0.0_dp
591 : fail = .TRUE.
592 : GOTO 300
593 : END IF
594 25 : step = (lam1 + lam2)/2
595 25 : fun3 = 0.0_dp
596 775 : DO indf = 1, ndf
597 775 : fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
598 : END DO
599 25 : fun3 = fun3 - step
600 :
601 25 : IF (ABS(step - lam2) < stol) THEN
602 : ln = step
603 : GOTO 200
604 : END IF
605 :
606 24 : IF (fun3*fun1 < stol) THEN
607 : lam2 = step
608 : ELSE
609 24 : lam1 = step
610 : END IF
611 : END DO
612 : END IF
613 : END DO
614 :
615 : 200 CONTINUE
616 802 : IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
617 : (eigval(iref) > 0.0_dp))) THEN
618 :
619 1 : IF (.NOT. bisec) GOTO 100
620 : ln = 0.0_dp
621 : lp = 0.0_dp
622 : fail = .TRUE.
623 : END IF
624 :
625 : 300 CONTINUE
626 :
627 801 : IF (fail .AND. .NOT. set) THEN
628 0 : set = .TRUE.
629 0 : DO indf = 1, ndf
630 0 : eigval(indf) = eigval(indf)*work(indf)
631 : END DO
632 : CYCLE
633 : END IF
634 :
635 801 : IF (.NOT. set) THEN
636 80871 : work(1:ndf) = one
637 : END IF
638 :
639 80871 : DO indf = 1, ndf
640 80871 : eigval(indf) = eigval(indf) - ln
641 : END DO
642 : EXIT
643 : END DO
644 :
645 801 : CALL timestop(handle)
646 :
647 801 : END SUBROUTINE rat_fun_opt
648 :
649 : ! **************************************************************************************************
650 : !> \brief ...
651 : !> \param ndf ...
652 : !> \param dx ...
653 : !> \param dg ...
654 : !> \param hess_mat ...
655 : !> \param work ...
656 : !> \param para_env ...
657 : ! **************************************************************************************************
658 6524 : SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env)
659 : INTEGER, INTENT(IN) :: ndf
660 : REAL(KIND=dp), INTENT(INOUT) :: dx(ndf), dg(ndf)
661 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
662 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
663 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
664 :
665 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bfgs'
666 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
667 :
668 : INTEGER :: handle, i, j, k, l, ncol_local, &
669 : nrow_local
670 3262 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
671 : REAL(KIND=dp) :: DDOT, dxw, gdx
672 3262 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_hes
673 :
674 3262 : CALL timeset(routineN, handle)
675 :
676 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
677 3262 : local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
678 :
679 254887 : work = zero
680 153028 : DO i = 1, nrow_local
681 149766 : j = row_indices(i)
682 23299840 : DO k = 1, ncol_local
683 23146812 : l = col_indices(k)
684 23296578 : work(j) = work(j) + local_hes(i, k)*dx(l)
685 : END DO
686 : END DO
687 :
688 3262 : CALL para_env%sum(work)
689 :
690 3262 : gdx = DDOT(ndf, dg, 1, dx, 1)
691 3262 : gdx = one/gdx
692 3262 : dxw = DDOT(ndf, dx, 1, work, 1)
693 3262 : dxw = one/dxw
694 :
695 153028 : DO i = 1, nrow_local
696 149766 : j = row_indices(i)
697 23299840 : DO k = 1, ncol_local
698 23146812 : l = col_indices(k)
699 : local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
700 23296578 : dxw*work(j)*work(l)
701 : END DO
702 : END DO
703 :
704 3262 : CALL timestop(handle)
705 :
706 3262 : END SUBROUTINE bfgs
707 :
708 : ! **************************************************************************************************
709 : !> \brief ...
710 : !> \param ndf ...
711 : !> \param eigval ...
712 : !> \param work ...
713 : ! **************************************************************************************************
714 801 : SUBROUTINE set_hes_eig(ndf, eigval, work)
715 : INTEGER, INTENT(IN) :: ndf
716 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf), work(ndf)
717 :
718 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_hes_eig'
719 : REAL(KIND=dp), PARAMETER :: max_neg = -0.5_dp, max_pos = 5.0_dp, &
720 : min_eig = 0.005_dp, one = 1.0_dp
721 :
722 : INTEGER :: handle, indf
723 : LOGICAL :: neg
724 :
725 801 : CALL timeset(routineN, handle)
726 :
727 80871 : DO indf = 1, ndf
728 80070 : IF (eigval(indf) < 0.0_dp) neg = .TRUE.
729 80871 : IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
730 : END DO
731 80871 : DO indf = 1, ndf
732 80871 : IF (eigval(indf) < 0.0_dp) THEN
733 1 : IF (eigval(indf) < max_neg) THEN
734 0 : eigval(indf) = max_neg
735 1 : ELSE IF (eigval(indf) > -min_eig) THEN
736 1 : eigval(indf) = -min_eig
737 : END IF
738 80069 : ELSE IF (eigval(indf) < 1000.0_dp) THEN
739 80069 : IF (eigval(indf) < min_eig) THEN
740 289 : eigval(indf) = min_eig
741 79780 : ELSE IF (eigval(indf) > max_pos) THEN
742 0 : eigval(indf) = max_pos
743 : END IF
744 : END IF
745 : END DO
746 :
747 80871 : DO indf = 1, ndf
748 80871 : IF (eigval(indf) < 0.0_dp) THEN
749 1 : work(indf) = -one
750 : ELSE
751 80069 : work(indf) = one
752 : END IF
753 : END DO
754 :
755 801 : CALL timestop(handle)
756 :
757 801 : END SUBROUTINE set_hes_eig
758 :
759 : ! **************************************************************************************************
760 : !> \brief ...
761 : !> \param ndf ...
762 : !> \param eigval ...
763 : !> \param eigvec_mat ...
764 : !> \param hess_tmp ...
765 : !> \param dr ...
766 : !> \param g ...
767 : !> \param para_env ...
768 : !> \param use_rfo ...
769 : ! **************************************************************************************************
770 12384 : SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
771 :
772 : INTEGER, INTENT(IN) :: ndf
773 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf)
774 : TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat, hess_tmp
775 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
776 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
777 : LOGICAL :: use_rfo
778 :
779 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
780 :
781 : INTEGER :: i, indf, j, k, l, ncol_local, nrow_local
782 4128 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
783 4128 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
784 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
785 : TYPE(cp_fm_type) :: tmp
786 :
787 4128 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
788 4128 : IF (use_rfo) THEN
789 80871 : DO indf = 1, ndf
790 80871 : eigval(indf) = one/eigval(indf)
791 : END DO
792 : ELSE
793 197826 : DO indf = 1, ndf
794 197826 : eigval(indf) = one/MAX(0.0001_dp, eigval(indf))
795 : END DO
796 : END IF
797 :
798 4128 : CALL cp_fm_column_scale(hess_tmp, eigval)
799 4128 : CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
800 4128 : CALL cp_fm_create(tmp, matrix_struct, name="tmp")
801 4128 : CALL cp_fm_set_all(tmp, alpha=zero)
802 :
803 4128 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
804 :
805 4128 : CALL cp_fm_transpose(tmp, hess_tmp)
806 4128 : CALL cp_fm_release(tmp)
807 :
808 : ! New step
809 :
810 : CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
811 4128 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
812 :
813 278697 : dr = 0.0_dp
814 167316 : DO i = 1, nrow_local
815 163188 : j = row_indices(i)
816 25014102 : DO k = 1, ncol_local
817 24846786 : l = col_indices(k)
818 25009974 : dr(j) = dr(j) - local_data(i, k)*g(l)
819 : END DO
820 : END DO
821 :
822 4128 : CALL para_env%sum(dr)
823 :
824 4128 : END SUBROUTINE geoopt_get_step
825 :
826 : ! **************************************************************************************************
827 : !> \brief ...
828 : !> \param ndf ...
829 : !> \param step ...
830 : !> \param rad ...
831 : !> \param rat ...
832 : !> \param dr ...
833 : !> \param output_unit ...
834 : ! **************************************************************************************************
835 4128 : SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
836 : INTEGER, INTENT(IN) :: ndf
837 : REAL(KIND=dp), INTENT(INOUT) :: step, rad, rat, dr(ndf)
838 : INTEGER, INTENT(IN) :: output_unit
839 :
840 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trust_radius'
841 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
842 :
843 : INTEGER :: handle
844 : REAL(KIND=dp) :: scal
845 :
846 4128 : CALL timeset(routineN, handle)
847 :
848 278697 : step = MAXVAL(ABS(dr))
849 4128 : scal = MAX(one, rad/step)
850 :
851 4128 : IF (step > rad) THEN
852 372 : rat = rad/step
853 372 : CALL DSCAL(ndf, rat, dr, 1)
854 372 : step = rad
855 372 : IF (output_unit > 0) THEN
856 : WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") &
857 189 : " Step is scaled; Scaling factor = ", rat
858 189 : CALL m_flush(output_unit)
859 : END IF
860 : END IF
861 4128 : CALL timestop(handle)
862 :
863 4128 : END SUBROUTINE trust_radius
864 :
865 : ! **************************************************************************************************
866 : !> \brief ...
867 : !> \param ndf ...
868 : !> \param work ...
869 : !> \param hess_mat ...
870 : !> \param dr ...
871 : !> \param g ...
872 : !> \param conv ...
873 : !> \param pred ...
874 : !> \param para_env ...
875 : ! **************************************************************************************************
876 8256 : SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
877 :
878 : INTEGER, INTENT(IN) :: ndf
879 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
880 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
881 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
882 : LOGICAL, INTENT(INOUT) :: conv
883 : REAL(KIND=dp), INTENT(INOUT) :: pred
884 : TYPE(mp_para_env_type), POINTER :: para_env
885 :
886 : CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_predict'
887 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
888 :
889 : INTEGER :: handle, i, j, k, l, ncol_local, &
890 : nrow_local
891 4128 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
892 : REAL(KIND=dp) :: DDOT, ener1, ener2
893 4128 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
894 :
895 4128 : CALL timeset(routineN, handle)
896 :
897 4128 : ener1 = DDOT(ndf, g, 1, dr, 1)
898 :
899 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
900 4128 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
901 :
902 278697 : work = zero
903 167316 : DO i = 1, nrow_local
904 163188 : j = row_indices(i)
905 25014102 : DO k = 1, ncol_local
906 24846786 : l = col_indices(k)
907 25009974 : work(j) = work(j) + local_data(i, k)*dr(l)
908 : END DO
909 : END DO
910 :
911 4128 : CALL para_env%sum(work)
912 4128 : ener2 = DDOT(ndf, dr, 1, work, 1)
913 4128 : pred = ener1 + 0.5_dp*ener2
914 4128 : conv = .FALSE.
915 4128 : CALL timestop(handle)
916 :
917 4128 : END SUBROUTINE energy_predict
918 :
919 : ! **************************************************************************************************
920 : !> \brief ...
921 : !> \param rat ...
922 : !> \param rad ...
923 : !> \param step ...
924 : !> \param ediff ...
925 : ! **************************************************************************************************
926 763 : SUBROUTINE update_trust_rad(rat, rad, step, ediff)
927 :
928 : REAL(KIND=dp), INTENT(INOUT) :: rat, rad, step, ediff
929 :
930 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_trust_rad'
931 : REAL(KIND=dp), PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp
932 :
933 : INTEGER :: handle
934 :
935 763 : CALL timeset(routineN, handle)
936 :
937 763 : IF (rat > 4.0_dp) THEN
938 0 : IF (ediff < 0.0_dp) THEN
939 0 : rad = step*0.5_dp
940 : ELSE
941 0 : rad = step*0.25_dp
942 : END IF
943 763 : ELSE IF (rat > 2.0_dp) THEN
944 0 : IF (ediff < 0.0_dp) THEN
945 0 : rad = step*0.75_dp
946 : ELSE
947 0 : rad = step*0.5_dp
948 : END IF
949 763 : ELSE IF (rat > 4.0_dp/3.0_dp) THEN
950 0 : IF (ediff < 0.0_dp) THEN
951 0 : rad = step
952 : ELSE
953 0 : rad = step*0.75_dp
954 : END IF
955 763 : ELSE IF (rat > 10.0_dp/9.0_dp) THEN
956 0 : IF (ediff < 0.0_dp) THEN
957 0 : rad = step*1.25_dp
958 : ELSE
959 0 : rad = step
960 : END IF
961 763 : ELSE IF (rat > 0.9_dp) THEN
962 73 : IF (ediff < 0.0_dp) THEN
963 73 : rad = step*1.5_dp
964 : ELSE
965 0 : rad = step*1.25_dp
966 : END IF
967 690 : ELSE IF (rat > 0.75_dp) THEN
968 97 : IF (ediff < 0.0_dp) THEN
969 94 : rad = step*1.25_dp
970 : ELSE
971 3 : rad = step
972 : END IF
973 593 : ELSE IF (rat > 0.5_dp) THEN
974 86 : IF (ediff < 0.0_dp) THEN
975 85 : rad = step
976 : ELSE
977 1 : rad = step*0.75_dp
978 : END IF
979 507 : ELSE IF (rat > 0.25_dp) THEN
980 5 : IF (ediff < 0.0_dp) THEN
981 5 : rad = step*0.75_dp
982 : ELSE
983 0 : rad = step*0.5_dp
984 : END IF
985 502 : ELSE IF (ediff < 0.0_dp) THEN
986 501 : rad = step*0.5_dp
987 : ELSE
988 1 : rad = step*0.25_dp
989 : END IF
990 :
991 763 : rad = MAX(rad, min_trust)
992 763 : rad = MIN(rad, max_trust)
993 763 : CALL timestop(handle)
994 :
995 763 : END SUBROUTINE update_trust_rad
996 :
997 : ! **************************************************************************************************
998 :
999 : ! **************************************************************************************************
1000 : !> \brief ...
1001 : !> \param geo_section ...
1002 : !> \param hess_mat ...
1003 : !> \param logger ...
1004 : ! **************************************************************************************************
1005 3539 : SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger)
1006 :
1007 : TYPE(section_vals_type), POINTER :: geo_section
1008 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1009 : TYPE(cp_logger_type), POINTER :: logger
1010 :
1011 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bfgs_hessian'
1012 :
1013 : INTEGER :: handle, hesunit
1014 :
1015 3539 : CALL timeset(routineN, handle)
1016 :
1017 : hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", &
1018 : extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", &
1019 3539 : file_position="REWIND")
1020 :
1021 3539 : CALL cp_fm_write_unformatted(hess_mat, hesunit)
1022 :
1023 3539 : CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
1024 :
1025 3539 : CALL timestop(handle)
1026 :
1027 3539 : END SUBROUTINE write_bfgs_hessian
1028 :
1029 : ! **************************************************************************************************
1030 : !> \brief Constructs model Hessian as described in https://doi.org/10.1016/0009-2614(95)00646-L.
1031 : !> \param force_env ...
1032 : !> \param hess_mat ...
1033 : !> \author Florian Schiffmann
1034 : ! **************************************************************************************************
1035 754 : SUBROUTINE construct_initial_hess(force_env, hess_mat)
1036 :
1037 : TYPE(force_env_type), POINTER :: force_env
1038 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1039 :
1040 : INTEGER :: i, iat_col, iat_row, iglobal, iind, j, &
1041 : jat_row, jglobal, jind, k, natom, &
1042 : ncol_local, nrow_local, z
1043 754 : INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row
1044 754 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1045 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij
1046 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij
1047 : REAL(KIND=dp), DIMENSION(3, 3) :: alpha, r0
1048 754 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: fixed, local_data
1049 : TYPE(cell_type), POINTER :: cell
1050 : TYPE(cp_subsys_type), POINTER :: subsys
1051 : TYPE(particle_list_type), POINTER :: particles
1052 :
1053 754 : CALL cite_reference(Lindh1995)
1054 :
1055 754 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1056 : CALL cp_subsys_get(subsys, &
1057 754 : particles=particles)
1058 :
1059 3016 : alpha(1, :) = [1._dp, 0.3949_dp, 0.3949_dp]
1060 3016 : alpha(2, :) = [0.3494_dp, 0.2800_dp, 0.2800_dp]
1061 3016 : alpha(3, :) = [0.3494_dp, 0.2800_dp, 0.1800_dp]
1062 :
1063 3016 : r0(1, :) = [1.35_dp, 2.10_dp, 2.53_dp]
1064 3016 : r0(2, :) = [2.10_dp, 2.87_dp, 3.40_dp]
1065 3016 : r0(3, :) = [2.53_dp, 3.40_dp, 3.40_dp]
1066 :
1067 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1068 754 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1069 754 : natom = particles%n_els
1070 2262 : ALLOCATE (at_row(natom))
1071 3016 : ALLOCATE (rho_ij(natom, natom))
1072 2262 : ALLOCATE (d_ij(natom, natom))
1073 3770 : ALLOCATE (r_ij(natom, natom, 3))
1074 2262 : ALLOCATE (fixed(3, natom))
1075 21642 : fixed = 1.0_dp
1076 754 : CALL fix_atom_control(force_env, fixed)
1077 3016 : DO i = 1, 3
1078 34348 : CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1079 : END DO
1080 754 : rho_ij = 0
1081 : !XXXX insert proper rows !XXX
1082 5976 : at_row = 3
1083 5976 : DO i = 1, natom
1084 5222 : CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
1085 5222 : IF (z <= 10) at_row(i) = 2
1086 11198 : IF (z <= 2) at_row(i) = 1
1087 : END DO
1088 5222 : DO i = 2, natom
1089 4468 : iat_row = at_row(i)
1090 64304 : DO j = 1, i - 1
1091 59082 : jat_row = at_row(j)
1092 : !pbc for a distance vector
1093 236328 : r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
1094 236328 : r_ij(i, j, :) = -r_ij(j, i, :)
1095 236328 : d_ij(j, i) = NORM2(r_ij(j, i, :))
1096 59082 : d_ij(i, j) = d_ij(j, i)
1097 59082 : rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1098 63550 : rho_ij(i, j) = rho_ij(j, i)
1099 : END DO
1100 : END DO
1101 16420 : DO i = 1, ncol_local
1102 15666 : iglobal = col_indices(i)
1103 15666 : iind = MOD(iglobal - 1, 3) + 1
1104 15666 : iat_col = (iglobal + 2)/3
1105 15666 : IF (iat_col > natom) CYCLE
1106 662287 : DO j = 1, nrow_local
1107 645867 : jglobal = row_indices(j)
1108 645867 : jind = MOD(jglobal - 1, 3) + 1
1109 645867 : iat_row = (jglobal + 2)/3
1110 645867 : IF (iat_row > natom) CYCLE
1111 645867 : IF (iat_row /= iat_col) THEN
1112 616518 : IF (d_ij(iat_row, iat_col) < 6.0_dp) THEN
1113 : local_data(j, i) = local_data(j, i) + &
1114 217422 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1115 : END IF
1116 : ELSE
1117 : local_data(j, i) = local_data(j, i) + &
1118 29349 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1119 : END IF
1120 645867 : IF (iat_col /= iat_row) THEN
1121 616518 : IF (d_ij(iat_row, iat_col) < 6.0_dp) THEN
1122 : local_data(j, i) = local_data(j, i) - &
1123 : dist_second_deriv(r_ij(iat_col, iat_row, :), &
1124 1521954 : iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1125 : END IF
1126 : ELSE
1127 675216 : DO k = 1, natom
1128 645867 : IF (k == iat_col) CYCLE
1129 645867 : IF (d_ij(iat_row, k) < 6.0_dp) THEN
1130 : local_data(j, i) = local_data(j, i) + &
1131 : dist_second_deriv(r_ij(iat_col, k, :), &
1132 1521954 : iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1133 : END IF
1134 : END DO
1135 : END IF
1136 661533 : IF (fixed(jind, iat_row) < 0.5_dp .OR. fixed(iind, iat_col) < 0.5_dp) THEN
1137 10161 : local_data(j, i) = 0.0_dp
1138 10161 : IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
1139 : END IF
1140 : END DO
1141 : END DO
1142 754 : DEALLOCATE (fixed)
1143 754 : DEALLOCATE (rho_ij)
1144 754 : DEALLOCATE (d_ij)
1145 754 : DEALLOCATE (r_ij)
1146 754 : DEALLOCATE (at_row)
1147 :
1148 1508 : END SUBROUTINE construct_initial_hess
1149 :
1150 : ! **************************************************************************************************
1151 : !> \brief ...
1152 : !> \param r1 ...
1153 : !> \param i ...
1154 : !> \param j ...
1155 : !> \param d ...
1156 : !> \param rho ...
1157 : !> \return ...
1158 : ! **************************************************************************************************
1159 434844 : FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv)
1160 : REAL(KIND=dp), DIMENSION(3) :: r1
1161 : INTEGER :: i, j
1162 : REAL(KIND=dp) :: d, rho, deriv
1163 :
1164 434844 : deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
1165 434844 : END FUNCTION dist_second_deriv
1166 :
1167 : ! **************************************************************************************************
1168 : !> \brief ...
1169 : !> \param r_ij ...
1170 : !> \param d_ij ...
1171 : !> \param rho_ij ...
1172 : !> \param idir ...
1173 : !> \param jdir ...
1174 : !> \param iat_der ...
1175 : !> \param jat_der ...
1176 : !> \param natom ...
1177 : !> \return ...
1178 : ! **************************************************************************************************
1179 246771 : FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv)
1180 : REAL(KIND=dp), DIMENSION(:, :, :) :: r_ij
1181 : REAL(KIND=dp), DIMENSION(:, :) :: d_ij, rho_ij
1182 : INTEGER :: idir, jdir, iat_der, jat_der, natom
1183 : REAL(KIND=dp) :: deriv
1184 :
1185 : INTEGER :: i, iat, idr, j, jat, jdr
1186 : REAL(KIND=dp) :: d12, d23, d31, D_mat(3, 2), denom1, &
1187 : denom2, denom3, ka1, ka2, ka3, rho12, &
1188 : rho23, rho31, rsst1, rsst2, rsst3
1189 : REAL(KIND=dp), DIMENSION(3) :: r12, r23, r31
1190 :
1191 246771 : deriv = 0._dp
1192 246771 : IF (iat_der == jat_der) THEN
1193 645867 : DO i = 1, natom - 1
1194 616518 : IF (rho_ij(iat_der, i) < 0.00001) CYCLE
1195 3408813 : DO j = i + 1, natom
1196 3183588 : IF (rho_ij(iat_der, j) < 0.00001) CYCLE
1197 947088 : IF (i == iat_der .OR. j == iat_der) CYCLE
1198 947088 : IF (iat_der < i .OR. iat_der > j) THEN
1199 5773320 : r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1200 577332 : d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1201 577332 : rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1202 : ELSE
1203 3697560 : r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1204 369756 : d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1205 369756 : rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1206 : END IF
1207 947088 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1208 9470880 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1209 947088 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1210 947088 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1211 947088 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1212 947088 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1213 947088 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1214 947088 : D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1215 947088 : D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1216 947088 : D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1217 947088 : D_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
1218 : D_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
1219 947088 : rsst3*r12(idir)/(d31*d12**3)
1220 : D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1221 947088 : rsst3*r12(jdir)/(d31*d12**3)
1222 947088 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1223 947088 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1224 947088 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1225 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1226 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1227 3800106 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1228 :
1229 : END DO
1230 : END DO
1231 : ELSE
1232 6123942 : DO i = 1, natom
1233 5906520 : IF (i == iat_der .OR. i == jat_der) CYCLE
1234 5471676 : IF (jat_der < iat_der) THEN
1235 2735838 : iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1236 : ELSE
1237 2735838 : iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1238 : END IF
1239 5471676 : IF (jat < i .OR. iat > i) THEN
1240 43236540 : r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1241 4323654 : d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1242 4323654 : rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1243 : ELSE
1244 11480220 : r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1245 1148022 : d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1246 1148022 : rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1247 : END IF
1248 5471676 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1249 54716760 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1250 5471676 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1251 5471676 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1252 5471676 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1253 5471676 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1254 5471676 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1255 5471676 : D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1256 5471676 : D_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
1257 : D_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
1258 5471676 : rsst3*r12(idr)/(d31*d12**3)
1259 5471676 : IF (jat < i .OR. iat > i) THEN
1260 : D_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
1261 4323654 : rsst1*r23(jdr)/(d12*d23**3)
1262 4323654 : D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1263 4323654 : D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1264 : ELSE
1265 1148022 : D_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
1266 : D_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
1267 1148022 : rsst2*r31(jdr)/(d23*d31**3)
1268 1148022 : D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1269 : END IF
1270 5471676 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1271 5471676 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1272 5471676 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1273 :
1274 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1275 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1276 6123942 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1277 : END DO
1278 : END IF
1279 246771 : deriv = 0.25_dp*deriv
1280 :
1281 246771 : END FUNCTION angle_second_deriv
1282 :
1283 : END MODULE bfgs_optimizer
|
