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