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