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