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 1115 : 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 1115 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
127 1115 : 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 1115 : NULLIFY (logger, g, blacs_env, spgr)
139 2230 : logger => cp_get_default_logger()
140 1115 : para_env => force_env%para_env
141 1115 : root_section => force_env%root_section
142 1115 : spgr => gopt_env%spgr
143 1115 : t_old = m_walltime()
144 :
145 1115 : CALL timeset(routineN, handle)
146 1115 : CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
147 1115 : print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
148 1115 : ionode = para_env%is_source()
149 1115 : maxiter = gopt_param%max_iter
150 1115 : conv = .FALSE.
151 1115 : rat = 0.0_dp
152 1115 : wildcard = " BFGS"
153 :
154 : ! Stop if not yet implemented
155 1115 : SELECT CASE (gopt_env%type_id)
156 : CASE (default_ts_method_id)
157 1115 : CPABORT("BFGS method not yet working with DIMER")
158 : END SELECT
159 :
160 1115 : CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
161 1115 : CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
162 1115 : 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 1115 : extension=".geoLog")
165 1115 : IF (output_unit > 0) THEN
166 575 : IF (use_rfo) THEN
167 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
168 0 : "BFGS| Use rational function optimization for step estimation: ", "YES"
169 : ELSE
170 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A3)") &
171 575 : "BFGS| Use rational function optimization for step estimation: ", " NO"
172 : END IF
173 575 : IF (use_mod_hes) THEN
174 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
175 507 : "BFGS| Use model Hessian for initial guess: ", "YES"
176 : ELSE
177 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
178 68 : "BFGS| Use model Hessian for initial guess: ", " NO"
179 : END IF
180 575 : 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 574 : "BFGS| Restart Hessian: ", " NO"
186 : END IF
187 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
188 575 : "BFGS| Trust radius: ", rad
189 : END IF
190 :
191 1115 : ndf = SIZE(x0)
192 1115 : nfree = gopt_env%nfree
193 1115 : 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 1115 : globenv%blacs_repeatable)
204 : CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
205 1115 : nrow_global=ndf, ncol_global=ndf)
206 1115 : CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
207 1115 : CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
208 1115 : CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
209 3345 : ALLOCATE (eigval(ndf))
210 62327 : eigval(:) = 0.0_dp
211 :
212 1115 : CALL force_env_get(force_env=force_env, subsys=subsys)
213 1115 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
214 1115 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
215 1115 : IF (use_mod_hes) THEN
216 979 : 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 979 : 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 1115 : IF (use_mod_hes) THEN
231 975 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=0.00_dp)
232 975 : CALL construct_initial_hess(gopt_env%force_env, hess_mat)
233 975 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
234 975 : 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 975 : 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 54453 : DO its = 1, SIZE(eigval)
242 54453 : IF (eigval(its) .LT. 0.1_dp) eigval(its) = 0.1_dp
243 : END DO
244 975 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
245 975 : CALL cp_fm_column_scale(eigvec_mat, eigval)
246 975 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
247 : END IF
248 : ELSE
249 140 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
250 : END IF
251 :
252 3345 : ALLOCATE (xold(ndf))
253 62327 : xold(:) = x0(:)
254 :
255 3345 : ALLOCATE (g(ndf))
256 62327 : g(:) = 0.0_dp
257 :
258 3345 : ALLOCATE (gold(ndf))
259 62327 : gold(:) = 0.0_dp
260 :
261 3345 : ALLOCATE (dx(ndf))
262 62327 : dx(:) = 0.0_dp
263 :
264 3345 : ALLOCATE (dg(ndf))
265 62327 : dg(:) = 0.0_dp
266 :
267 3345 : ALLOCATE (work(ndf))
268 62327 : work(:) = 0.0_dp
269 :
270 3345 : ALLOCATE (dr(ndf))
271 62327 : dr(:) = 0.0_dp
272 :
273 : ! find space_group
274 1115 : CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
275 1115 : 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 1115 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
283 1115 : 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 1115 : .FALSE., gopt_env%force_env%para_env)
288 :
289 : ! Symmetrize coordinates and forces
290 1115 : 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 1115 : emin = etot
297 1115 : t_now = m_walltime()
298 1115 : t_diff = t_now - t_old
299 1115 : t_old = t_now
300 1115 : CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
301 6194 : DO its = iter_nr + 1, maxiter
302 6184 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
303 6184 : CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
304 6184 : CALL gopt_f_ii(its, output_unit)
305 :
306 : ! Hessian update/restarting
307 6184 : 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 6182 : IF ((its - iter_nr) > 1) THEN
317 : ! Symmetrize old coordinates and old forces
318 5079 : 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 519438 : DO indf = 1, ndf
324 514359 : dx(indf) = x0(indf) - xold(indf)
325 519438 : dg(indf) = g(indf) - gold(indf)
326 : END DO
327 :
328 5079 : CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
329 :
330 : ! Symmetrize coordinates and forces change
331 5079 : 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 5079 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
338 4091 : 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 6184 : 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 581611 : xold(:) = x0
351 581611 : gold(:) = g
352 :
353 : ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
354 6184 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
355 :
356 6184 : 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 6184 : 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 6184 : IF (use_rfo) THEN
368 617 : CALL set_hes_eig(ndf, eigval, work)
369 19127 : dx(:) = eigval
370 617 : CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
371 : END IF
372 6184 : CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
373 :
374 : ! Symmetrize dr
375 6184 : IF (spgr%keep_space_group) THEN
376 4 : CALL spgr_apply_rotations_force(spgr, dr)
377 : END IF
378 :
379 6184 : CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
380 :
381 : ! Update the atomic positions
382 581611 : x0 = x0 + dr
383 :
384 : ! Symmetrize coordinates
385 6184 : IF (spgr%keep_space_group) THEN
386 4 : CALL spgr_apply_rotations_coord(spgr, x0)
387 : END IF
388 :
389 6184 : CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
390 6184 : 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 6184 : .FALSE., gopt_env%force_env%para_env)
395 :
396 6184 : ediff = etot - eold
397 :
398 : ! Symmetrize forces
399 6184 : 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 6184 : CALL external_control(should_stop, "GEO", globenv=globenv)
405 6184 : IF (should_stop) EXIT
406 :
407 : ! Some IO and Convergence check
408 6184 : t_now = m_walltime()
409 6184 : t_diff = t_now - t_old
410 6184 : 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 6184 : step, rad, used_time=t_diff)
414 :
415 6184 : IF (conv .OR. (its == maxiter)) EXIT
416 5079 : IF (etot < emin) emin = etot
417 18562 : IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
418 : END DO
419 :
420 1115 : IF (its == maxiter .AND. (.NOT. conv)) THEN
421 561 : CALL print_geo_opt_nc(gopt_env, output_unit)
422 : END IF
423 :
424 : ! Write final information, if converged
425 1115 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
426 1115 : 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 1115 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
429 :
430 1115 : CALL cp_fm_struct_release(fm_struct_hes)
431 1115 : CALL cp_fm_release(hess_mat)
432 1115 : CALL cp_fm_release(eigvec_mat)
433 1115 : CALL cp_fm_release(hess_tmp)
434 :
435 1115 : CALL cp_blacs_env_release(blacs_env)
436 1115 : DEALLOCATE (xold)
437 1115 : DEALLOCATE (g)
438 1115 : DEALLOCATE (gold)
439 1115 : DEALLOCATE (dx)
440 1115 : DEALLOCATE (dg)
441 1115 : DEALLOCATE (eigval)
442 1115 : DEALLOCATE (work)
443 1115 : DEALLOCATE (dr)
444 :
445 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
446 1115 : "PRINT%PROGRAM_RUN_INFO")
447 1115 : CALL timestop(handle)
448 :
449 4460 : 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 1234 : 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 617 : 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 617 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
479 :
480 617 : CALL timeset(routineN, handle)
481 :
482 617 : stol = 1.0E-8_dp
483 617 : ssize = 0.2_dp
484 617 : maxit = 999
485 617 : fail = .FALSE.
486 617 : bisec = .FALSE.
487 :
488 19127 : dg = 0._dp
489 :
490 : CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
491 617 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
492 :
493 19127 : DO i = 1, nrow_local
494 18510 : j = row_indices(i)
495 574427 : DO k = 1, ncol_local
496 555300 : l = col_indices(k)
497 573810 : dg(l) = dg(l) + local_data(i, k)*g(j)
498 : END DO
499 : END DO
500 617 : CALL para_env%sum(dg)
501 :
502 617 : set = .FALSE.
503 :
504 : DO
505 :
506 : ! calculating Lambda
507 :
508 617 : lp = 0.0_dp
509 617 : iref = 1
510 617 : ln = 0.0_dp
511 617 : IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
512 :
513 617 : iter = 0
514 : DO
515 1983 : iter = iter + 1
516 1983 : fun = 0.0_dp
517 1983 : fung = 0.0_dp
518 61473 : DO indf = 1, ndf
519 59490 : fun = fun + dg(indf)**2/(ln - eigval(indf))
520 61473 : fung = fung - dg(indf)**2/(ln - eigval(indf)**2)
521 : END DO
522 1983 : fun = fun - ln
523 1983 : fung = fung - one
524 1983 : step = fun/fung
525 1983 : ln = ln - step
526 1983 : IF (ABS(step) < stol) GOTO 200
527 1367 : IF (iter >= maxit) EXIT
528 : END DO
529 : 100 CONTINUE
530 1 : bisec = .TRUE.
531 1 : iter = 0
532 1 : maxit = 9999
533 1 : lam1 = 0.0_dp
534 1 : IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
535 : fun1 = 0.0_dp
536 31 : DO indf = 1, ndf
537 31 : fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
538 : END DO
539 1 : fun1 = fun1 - lam1
540 1 : step = ABS(lam1)/1000.0_dp
541 1 : IF (step < ssize) step = ssize
542 : DO
543 1 : iter = iter + 1
544 1 : IF (iter > maxit) THEN
545 : ln = 0.0_dp
546 617 : lp = 0.0_dp
547 : fail = .TRUE.
548 : GOTO 300
549 : END IF
550 1 : fun2 = 0.0_dp
551 1 : lam2 = lam1 - iter*step
552 31 : DO indf = 1, ndf
553 31 : fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf))
554 : END DO
555 1 : fun2 = fun2 - lam2
556 618 : IF (fun2*fun1 < 0.0_dp) THEN
557 : iter = 0
558 : DO
559 25 : iter = iter + 1
560 25 : IF (iter > maxit) THEN
561 : ln = 0.0_dp
562 : lp = 0.0_dp
563 : fail = .TRUE.
564 : GOTO 300
565 : END IF
566 25 : step = (lam1 + lam2)/2
567 25 : fun3 = 0.0_dp
568 775 : DO indf = 1, ndf
569 775 : fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
570 : END DO
571 25 : fun3 = fun3 - step
572 :
573 25 : IF (ABS(step - lam2) < stol) THEN
574 : ln = step
575 : GOTO 200
576 : END IF
577 :
578 24 : IF (fun3*fun1 < stol) THEN
579 : lam2 = step
580 : ELSE
581 24 : lam1 = step
582 : END IF
583 : END DO
584 : END IF
585 : END DO
586 :
587 : 200 CONTINUE
588 618 : IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
589 : (eigval(iref) > 0.0_dp))) THEN
590 :
591 1 : 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 617 : 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 617 : IF (.NOT. set) THEN
608 19127 : work(1:ndf) = one
609 : END IF
610 :
611 19127 : DO indf = 1, ndf
612 19127 : eigval(indf) = eigval(indf) - ln
613 : END DO
614 : EXIT
615 : END DO
616 :
617 617 : CALL timestop(handle)
618 :
619 617 : 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 10158 : 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 5079 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
643 : REAL(KIND=dp) :: DDOT, dxw, gdx
644 5079 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_hes
645 :
646 5079 : CALL timeset(routineN, handle)
647 :
648 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
649 5079 : local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
650 :
651 519438 : work = zero
652 286107 : DO i = 1, nrow_local
653 281028 : j = row_indices(i)
654 43170963 : DO k = 1, ncol_local
655 42884856 : l = col_indices(k)
656 43165884 : work(j) = work(j) + local_hes(i, k)*dx(l)
657 : END DO
658 : END DO
659 :
660 5079 : CALL para_env%sum(work)
661 :
662 5079 : gdx = DDOT(ndf, dg, 1, dx, 1)
663 5079 : gdx = one/gdx
664 5079 : dxw = DDOT(ndf, dx, 1, work, 1)
665 5079 : dxw = one/dxw
666 :
667 286107 : DO i = 1, nrow_local
668 281028 : j = row_indices(i)
669 43170963 : DO k = 1, ncol_local
670 42884856 : l = col_indices(k)
671 : local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
672 43165884 : dxw*work(j)*work(l)
673 : END DO
674 : END DO
675 :
676 5079 : CALL timestop(handle)
677 :
678 5079 : END SUBROUTINE bfgs
679 :
680 : ! **************************************************************************************************
681 : !> \brief ...
682 : !> \param ndf ...
683 : !> \param eigval ...
684 : !> \param work ...
685 : ! **************************************************************************************************
686 617 : 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 617 : CALL timeset(routineN, handle)
698 :
699 19127 : DO indf = 1, ndf
700 18510 : IF (eigval(indf) < 0.0_dp) neg = .TRUE.
701 19127 : IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
702 : END DO
703 19127 : DO indf = 1, ndf
704 19127 : IF (eigval(indf) < 0.0_dp) THEN
705 1 : IF (eigval(indf) < max_neg) THEN
706 0 : eigval(indf) = max_neg
707 1 : ELSE IF (eigval(indf) > -min_eig) THEN
708 1 : eigval(indf) = -min_eig
709 : END IF
710 18509 : ELSE IF (eigval(indf) < 1000.0_dp) THEN
711 18509 : IF (eigval(indf) < min_eig) THEN
712 315 : eigval(indf) = min_eig
713 18194 : ELSE IF (eigval(indf) > max_pos) THEN
714 0 : eigval(indf) = max_pos
715 : END IF
716 : END IF
717 : END DO
718 :
719 19127 : DO indf = 1, ndf
720 19127 : IF (eigval(indf) < 0.0_dp) THEN
721 1 : work(indf) = -one
722 : ELSE
723 18509 : work(indf) = one
724 : END IF
725 : END DO
726 :
727 617 : CALL timestop(handle)
728 :
729 617 : 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 12368 : 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 6184 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
755 6184 : 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 6184 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
760 6184 : IF (use_rfo) THEN
761 19127 : DO indf = 1, ndf
762 19127 : eigval(indf) = one/eigval(indf)
763 : END DO
764 : ELSE
765 562484 : DO indf = 1, ndf
766 562484 : eigval(indf) = one/MAX(0.0001_dp, eigval(indf))
767 : END DO
768 : END IF
769 :
770 6184 : CALL cp_fm_column_scale(hess_tmp, eigval)
771 6184 : CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
772 6184 : CALL cp_fm_create(tmp, matrix_struct, name="tmp")
773 :
774 6184 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
775 :
776 6184 : CALL cp_fm_transpose(tmp, hess_tmp)
777 6184 : 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 6184 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
783 :
784 581611 : dr = 0.0_dp
785 319651 : DO i = 1, nrow_local
786 313467 : j = row_indices(i)
787 48040918 : DO k = 1, ncol_local
788 47721267 : l = col_indices(k)
789 48034734 : dr(j) = dr(j) - local_data(i, k)*g(l)
790 : END DO
791 : END DO
792 :
793 6184 : CALL para_env%sum(dr)
794 :
795 6184 : 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 6184 : 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 6184 : CALL timeset(routineN, handle)
818 :
819 587795 : step = MAXVAL(ABS(dr))
820 6184 : scal = MAX(one, rad/step)
821 :
822 6184 : IF (step > rad) THEN
823 414 : rat = rad/step
824 414 : CALL DSCAL(ndf, rat, dr, 1)
825 414 : step = rad
826 414 : IF (output_unit > 0) THEN
827 : WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") &
828 211 : " Step is scaled; Scaling factor = ", rat
829 211 : CALL m_flush(output_unit)
830 : END IF
831 : END IF
832 6184 : CALL timestop(handle)
833 :
834 6184 : 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 12368 : 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 6184 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
863 : REAL(KIND=dp) :: DDOT, ener1, ener2
864 6184 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
865 :
866 6184 : CALL timeset(routineN, handle)
867 :
868 6184 : 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 6184 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
872 :
873 581611 : work = zero
874 319651 : DO i = 1, nrow_local
875 313467 : j = row_indices(i)
876 48040918 : DO k = 1, ncol_local
877 47721267 : l = col_indices(k)
878 48034734 : work(j) = work(j) + local_data(i, k)*dr(l)
879 : END DO
880 : END DO
881 :
882 6184 : CALL para_env%sum(work)
883 6184 : ener2 = DDOT(ndf, dr, 1, work, 1)
884 6184 : pred = ener1 + 0.5_dp*ener2
885 6184 : conv = .FALSE.
886 6184 : CALL timestop(handle)
887 :
888 6184 : END SUBROUTINE energy_predict
889 :
890 : ! **************************************************************************************************
891 : !> \brief ...
892 : !> \param rat ...
893 : !> \param rad ...
894 : !> \param step ...
895 : !> \param ediff ...
896 : ! **************************************************************************************************
897 592 : 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 592 : CALL timeset(routineN, handle)
907 :
908 592 : 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 592 : 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 592 : 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 592 : 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 592 : ELSE IF (rat > 0.9_dp) THEN
933 73 : IF (ediff < 0.0_dp) THEN
934 69 : rad = step*1.5_dp
935 : ELSE
936 4 : rad = step*1.25_dp
937 : END IF
938 519 : ELSE IF (rat > 0.75_dp) THEN
939 9 : IF (ediff < 0.0_dp) THEN
940 7 : rad = step*1.25_dp
941 : ELSE
942 2 : rad = step
943 : END IF
944 510 : ELSE IF (rat > 0.5_dp) THEN
945 75 : IF (ediff < 0.0_dp) THEN
946 75 : rad = step
947 : ELSE
948 0 : rad = step*0.75_dp
949 : END IF
950 435 : ELSE IF (rat > 0.25_dp) THEN
951 3 : IF (ediff < 0.0_dp) THEN
952 3 : rad = step*0.75_dp
953 : ELSE
954 0 : rad = step*0.5_dp
955 : END IF
956 432 : ELSE IF (ediff < 0.0_dp) THEN
957 430 : rad = step*0.5_dp
958 : ELSE
959 2 : rad = step*0.25_dp
960 : END IF
961 :
962 592 : rad = MAX(rad, min_trust)
963 592 : rad = MIN(rad, max_trust)
964 592 : CALL timestop(handle)
965 :
966 592 : END SUBROUTINE update_trust_rad
967 :
968 : ! **************************************************************************************************
969 :
970 : ! **************************************************************************************************
971 : !> \brief ...
972 : !> \param geo_section ...
973 : !> \param hess_mat ...
974 : !> \param logger ...
975 : ! **************************************************************************************************
976 5206 : 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 5206 : 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 5206 : file_position="REWIND")
991 :
992 5206 : CALL cp_fm_write_unformatted(hess_mat, hesunit)
993 :
994 5206 : CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
995 :
996 5206 : CALL timestop(handle)
997 :
998 5206 : END SUBROUTINE write_bfgs_hessian
999 : ! **************************************************************************************************
1000 :
1001 : ! **************************************************************************************************
1002 : !> \brief ...
1003 : !> \param force_env ...
1004 : !> \param hess_mat ...
1005 : ! **************************************************************************************************
1006 975 : 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 975 : INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row
1015 975 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1016 975 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij
1017 975 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij
1018 : REAL(KIND=dp), DIMENSION(3, 3) :: alpha, r0
1019 975 : 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 975 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1025 : CALL cp_subsys_get(subsys, &
1026 975 : particles=particles)
1027 :
1028 3900 : alpha(1, :) = (/1._dp, 0.3949_dp, 0.3949_dp/)
1029 3900 : alpha(2, :) = (/0.3494_dp, 0.2800_dp, 0.2800_dp/)
1030 3900 : alpha(3, :) = (/0.3494_dp, 0.2800_dp, 0.1800_dp/)
1031 :
1032 3900 : r0(1, :) = (/1.35_dp, 2.10_dp, 2.53_dp/)
1033 3900 : r0(2, :) = (/2.10_dp, 2.87_dp, 3.40_dp/)
1034 3900 : 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 975 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1038 975 : natom = particles%n_els
1039 2925 : ALLOCATE (at_row(natom))
1040 3900 : ALLOCATE (rho_ij(natom, natom))
1041 3900 : ALLOCATE (d_ij(natom, natom))
1042 4875 : ALLOCATE (r_ij(natom, natom, 3))
1043 2925 : ALLOCATE (fixed(3, natom))
1044 72279 : fixed = 1.0_dp
1045 975 : CALL fix_atom_control(force_env, fixed)
1046 3900 : DO i = 1, 3
1047 110856 : CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1048 : END DO
1049 916081 : rho_ij = 0
1050 : !XXXX insert proper rows !XXX
1051 18801 : at_row = 3
1052 18801 : DO i = 1, natom
1053 17826 : CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
1054 17826 : IF (z .LE. 10) at_row(i) = 2
1055 36627 : IF (z .LE. 2) at_row(i) = 1
1056 : END DO
1057 17826 : DO i = 2, natom
1058 16851 : iat_row = at_row(i)
1059 457553 : DO j = 1, i - 1
1060 439727 : jat_row = at_row(j)
1061 : !pbc for a distance vector
1062 1758908 : r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
1063 1758908 : r_ij(i, j, :) = -r_ij(j, i, :)
1064 1758908 : d_ij(j, i) = SQRT(DOT_PRODUCT(r_ij(j, i, :), r_ij(j, i, :)))
1065 439727 : d_ij(i, j) = d_ij(j, i)
1066 439727 : rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1067 456578 : rho_ij(i, j) = rho_ij(j, i)
1068 : END DO
1069 : END DO
1070 54453 : DO i = 1, ncol_local
1071 53478 : iglobal = col_indices(i)
1072 53478 : iind = MOD(iglobal - 1, 3) + 1
1073 53478 : iat_col = (iglobal + 2)/3
1074 53478 : IF (iat_col .GT. natom) CYCLE
1075 4181493 : DO j = 1, nrow_local
1076 4127040 : jglobal = row_indices(j)
1077 4127040 : jind = MOD(jglobal - 1, 3) + 1
1078 4127040 : iat_row = (jglobal + 2)/3
1079 4127040 : IF (iat_row .GT. natom) CYCLE
1080 4127040 : IF (iat_row .NE. iat_col) THEN
1081 4041108 : IF (d_ij(iat_row, iat_col) .LT. 6.0_dp) &
1082 : local_data(j, i) = local_data(j, i) + &
1083 343908 : 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 85932 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1087 : END IF
1088 4127040 : IF (iat_col .NE. iat_row) THEN
1089 4041108 : 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 2407356 : iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1093 : ELSE
1094 4212972 : DO k = 1, natom
1095 4127040 : IF (k == iat_col) CYCLE
1096 4041108 : 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 2493288 : iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1100 : END DO
1101 : END IF
1102 4180518 : 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 975 : DEALLOCATE (fixed)
1109 975 : DEALLOCATE (rho_ij)
1110 975 : DEALLOCATE (d_ij)
1111 975 : DEALLOCATE (r_ij)
1112 975 : DEALLOCATE (at_row)
1113 :
1114 2925 : 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 687816 : 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 687816 : deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
1131 687816 : 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 429840 : 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 429840 : deriv = 0._dp
1158 429840 : IF (iat_der == jat_der) THEN
1159 4127040 : DO i = 1, natom - 1
1160 4041108 : IF (rho_ij(iat_der, i) .LT. 0.00001) CYCLE
1161 25704207 : DO j = i + 1, natom
1162 24678504 : IF (rho_ij(iat_der, j) .LT. 0.00001) CYCLE
1163 6052374 : IF (i == iat_der .OR. j == iat_der) CYCLE
1164 6052374 : IF (iat_der .LT. i .OR. iat_der .GT. j) THEN
1165 38921670 : r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1166 3892167 : d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1167 3892167 : rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1168 : ELSE
1169 21602070 : r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1170 2160207 : d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1171 2160207 : rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1172 : END IF
1173 6052374 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1174 60523740 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1175 6052374 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1176 6052374 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1177 6052374 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1178 6052374 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1179 6052374 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1180 6052374 : D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1181 6052374 : D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1182 6052374 : D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1183 6052374 : 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 6052374 : rsst3*r12(idir)/(d31*d12**3)
1186 : D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1187 6052374 : rsst3*r12(jdir)/(d31*d12**3)
1188 6052374 : IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp
1189 6052374 : IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp
1190 6052374 : 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 28719612 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1194 :
1195 : END DO
1196 : END DO
1197 : ELSE
1198 10310040 : DO i = 1, natom
1199 9966132 : IF (i == iat_der .OR. i == jat_der) CYCLE
1200 9278316 : IF (jat_der .LT. iat_der) THEN
1201 4639158 : iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1202 : ELSE
1203 4639158 : iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1204 : END IF
1205 9278316 : IF (jat .LT. i .OR. iat .GT. i) THEN
1206 70107300 : r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1207 7010730 : d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1208 7010730 : rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1209 : ELSE
1210 22675860 : r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1211 2267586 : d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1212 2267586 : rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1213 : END IF
1214 9278316 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1215 92783160 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1216 9278316 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1217 9278316 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1218 9278316 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1219 9278316 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1220 9278316 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1221 9278316 : D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1222 9278316 : 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 9278316 : rsst3*r12(idr)/(d31*d12**3)
1225 9278316 : 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 7010730 : rsst1*r23(jdr)/(d12*d23**3)
1228 7010730 : D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1229 7010730 : D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1230 : ELSE
1231 2267586 : 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 2267586 : rsst2*r31(jdr)/(d23*d31**3)
1234 2267586 : D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1235 : END IF
1236 9278316 : IF (ABS(denom1) .LE. 0.011_dp) D_mat(1, 1) = 0.0_dp
1237 9278316 : IF (ABS(denom2) .LE. 0.011_dp) D_mat(2, 1) = 0.0_dp
1238 9278316 : 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 10310040 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1243 : END DO
1244 : END IF
1245 429840 : deriv = 0.25_dp*deriv
1246 :
1247 429840 : END FUNCTION angle_second_deriv
1248 :
1249 : END MODULE bfgs_optimizer
|
