Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 1239 : 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 1239 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
127 1239 : 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 1239 : NULLIFY (logger, g, blacs_env, spgr)
139 2478 : logger => cp_get_default_logger()
140 1239 : para_env => force_env%para_env
141 1239 : root_section => force_env%root_section
142 1239 : spgr => gopt_env%spgr
143 1239 : t_old = m_walltime()
144 :
145 1239 : CALL timeset(routineN, handle)
146 1239 : CALL section_vals_val_get(geo_section, "BFGS%TRUST_RADIUS", r_val=rad)
147 1239 : print_key => section_vals_get_subs_vals(geo_section, "BFGS%RESTART")
148 1239 : ionode = para_env%is_source()
149 1239 : maxiter = gopt_param%max_iter
150 1239 : conv = .FALSE.
151 1239 : rat = 0.0_dp
152 1239 : wildcard = " BFGS"
153 1239 : hes_filename = ""
154 :
155 : ! Stop if not yet implemented
156 1239 : SELECT CASE (gopt_env%type_id)
157 : CASE (default_ts_method_id)
158 1239 : CPABORT("BFGS method not yet working with DIMER")
159 : END SELECT
160 :
161 1239 : CALL section_vals_val_get(geo_section, "BFGS%USE_RAT_FUN_OPT", l_val=use_rfo)
162 1239 : CALL section_vals_val_get(geo_section, "BFGS%USE_MODEL_HESSIAN", l_val=use_mod_hes)
163 1239 : 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 1239 : extension=".geoLog")
166 1239 : IF (output_unit > 0) THEN
167 637 : 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 602 : "BFGS| Use rational function optimization for step estimation: ", " NO"
173 : END IF
174 637 : IF (use_mod_hes) THEN
175 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
176 531 : "BFGS| Use model Hessian for initial guess: ", "YES"
177 : ELSE
178 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
179 106 : "BFGS| Use model Hessian for initial guess: ", " NO"
180 : END IF
181 637 : 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 636 : "BFGS| Restart Hessian: ", " NO"
187 : END IF
188 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
189 637 : "BFGS| Trust radius: ", rad
190 : END IF
191 :
192 1239 : ndf = SIZE(x0)
193 1239 : nfree = gopt_env%nfree
194 1239 : 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 1239 : globenv%blacs_repeatable)
205 : CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
206 1239 : nrow_global=ndf, ncol_global=ndf)
207 1239 : CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
208 1239 : CALL cp_fm_create(hess_tmp, fm_struct_hes, name="hess_tmp")
209 1239 : CALL cp_fm_create(eigvec_mat, fm_struct_hes, name="eigvec_mat")
210 3717 : ALLOCATE (eigval(ndf))
211 88779 : eigval(:) = 0.0_dp
212 :
213 1239 : CALL force_env_get(force_env=force_env, subsys=subsys)
214 1239 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
215 1239 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, shell_present=shell_present)
216 1239 : IF (use_mod_hes) THEN
217 1027 : 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 1027 : 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 1239 : IF (use_mod_hes) THEN
232 1023 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=0.00_dp)
233 1023 : CALL construct_initial_hess(gopt_env%force_env, hess_mat)
234 1023 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
235 1023 : 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 1023 : 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 54837 : DO its = 1, SIZE(eigval)
243 54837 : IF (eigval(its) < 0.1_dp) eigval(its) = 0.1_dp
244 : END DO
245 1023 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
246 1023 : CALL cp_fm_column_scale(eigvec_mat, eigval)
247 1023 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
248 : END IF
249 : ELSE
250 216 : CALL cp_fm_set_all(hess_mat, alpha=zero, beta=one)
251 : END IF
252 :
253 2478 : ALLOCATE (xold(ndf))
254 88779 : xold(:) = x0(:)
255 :
256 2478 : ALLOCATE (g(ndf))
257 88779 : g(:) = 0.0_dp
258 :
259 2478 : ALLOCATE (gold(ndf))
260 88779 : gold(:) = 0.0_dp
261 :
262 2478 : ALLOCATE (dx(ndf))
263 88779 : dx(:) = 0.0_dp
264 :
265 2478 : ALLOCATE (dg(ndf))
266 88779 : dg(:) = 0.0_dp
267 :
268 2478 : ALLOCATE (work(ndf))
269 88779 : work(:) = 0.0_dp
270 :
271 2478 : ALLOCATE (dr(ndf))
272 88779 : dr(:) = 0.0_dp
273 :
274 : ! find space_group
275 1239 : CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
276 1239 : 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 1239 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
284 1239 : 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 1239 : .FALSE., gopt_env%force_env%para_env)
289 :
290 : ! Symmetrize coordinates and forces
291 1239 : 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 1239 : emin = etot
298 1239 : t_now = m_walltime()
299 1239 : t_diff = t_now - t_old
300 1239 : t_old = t_now
301 1239 : CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
302 6559 : DO its = iter_nr + 1, maxiter
303 6549 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
304 6549 : CALL section_vals_val_set(geo_section, "STEP_START_VAL", i_val=its)
305 6549 : CALL gopt_f_ii(its, output_unit)
306 :
307 : ! Hessian update/restarting
308 6549 : 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 6547 : IF ((its - iter_nr) > 1) THEN
330 : ! Symmetrize old coordinates and old forces
331 5320 : 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 582085 : DO indf = 1, ndf
337 576765 : dx(indf) = x0(indf) - xold(indf)
338 582085 : dg(indf) = g(indf) - gold(indf)
339 : END DO
340 :
341 5320 : CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
342 :
343 : ! Symmetrize coordinates and forces change
344 5320 : 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 5320 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
351 3781 : 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 6549 : 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 670710 : xold(:) = x0
364 670710 : gold(:) = g
365 :
366 : ! Copying hessian hes to (ndf x ndf) matrix hes_mat for diagonalization
367 6549 : CALL cp_fm_to_fm(hess_mat, hess_tmp)
368 :
369 6549 : 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 6549 : 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 6549 : IF (use_rfo) THEN
381 872 : CALL set_hes_eig(ndf, eigval, work)
382 105260 : dx(:) = eigval
383 872 : CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
384 : END IF
385 6549 : CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
386 :
387 : ! Symmetrize dr
388 6549 : IF (spgr%keep_space_group) THEN
389 4 : CALL spgr_apply_rotations_force(spgr, dr)
390 : END IF
391 :
392 6549 : CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
393 :
394 : ! Update the atomic positions
395 670710 : x0 = x0 + dr
396 :
397 : ! Symmetrize coordinates
398 6549 : IF (spgr%keep_space_group) THEN
399 4 : CALL spgr_apply_rotations_coord(spgr, x0)
400 : END IF
401 :
402 6549 : CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
403 6549 : 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 6549 : .FALSE., gopt_env%force_env%para_env)
408 :
409 6549 : ediff = etot - eold
410 :
411 : ! Symmetrize forces
412 6549 : 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 6549 : CALL external_control(should_stop, "GEO", globenv=globenv)
418 6549 : IF (should_stop) EXIT
419 :
420 : ! Some IO and Convergence check
421 6549 : t_now = m_walltime()
422 6549 : t_diff = t_now - t_old
423 6549 : 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 6549 : step, rad, used_time=t_diff)
427 :
428 6549 : IF (conv .OR. (its == maxiter)) EXIT
429 5320 : IF (etot < emin) emin = etot
430 19657 : IF (use_rfo) CALL update_trust_rad(rat, rad, step, ediff)
431 : END DO
432 :
433 1239 : IF (its == maxiter .AND. (.NOT. conv)) THEN
434 611 : CALL print_geo_opt_nc(gopt_env, output_unit)
435 : END IF
436 :
437 : ! show space_group
438 1239 : CALL section_vals_val_get(geo_section, "SHOW_SPACE_GROUP", l_val=spgr%show_space_group)
439 1239 : IF (spgr%show_space_group) THEN
440 2 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
441 2 : CALL print_spgr(spgr)
442 : END IF
443 :
444 : ! Write final information, if converged
445 1239 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
446 1239 : CALL write_bfgs_hessian(geo_section, hess_mat, logger)
447 : CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
448 1239 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
449 :
450 1239 : CALL cp_fm_struct_release(fm_struct_hes)
451 1239 : CALL cp_fm_release(hess_mat)
452 1239 : CALL cp_fm_release(eigvec_mat)
453 1239 : CALL cp_fm_release(hess_tmp)
454 :
455 1239 : CALL cp_blacs_env_release(blacs_env)
456 1239 : DEALLOCATE (xold)
457 1239 : DEALLOCATE (g)
458 1239 : DEALLOCATE (gold)
459 1239 : DEALLOCATE (dx)
460 1239 : DEALLOCATE (dg)
461 1239 : DEALLOCATE (eigval)
462 1239 : DEALLOCATE (work)
463 1239 : DEALLOCATE (dr)
464 :
465 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
466 1239 : "PRINT%PROGRAM_RUN_INFO")
467 1239 : CALL timestop(handle)
468 :
469 6195 : END SUBROUTINE geoopt_bfgs
470 :
471 : ! **************************************************************************************************
472 : !> \brief ...
473 : !> \param ndf ...
474 : !> \param dg ...
475 : !> \param eigval ...
476 : !> \param work ...
477 : !> \param eigvec_mat ...
478 : !> \param g ...
479 : !> \param para_env ...
480 : ! **************************************************************************************************
481 1744 : SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
482 :
483 : INTEGER, INTENT(IN) :: ndf
484 : REAL(KIND=dp), INTENT(INOUT) :: dg(ndf), eigval(ndf), work(ndf)
485 : TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat
486 : REAL(KIND=dp), INTENT(INOUT) :: g(ndf)
487 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
488 :
489 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rat_fun_opt'
490 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
491 :
492 : INTEGER :: handle, i, indf, iref, iter, j, k, l, &
493 : maxit, ncol_local, nrow_local
494 872 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
495 : LOGICAL :: bisec, fail, set
496 : REAL(KIND=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
497 : ln, lp, ssize, step, stol
498 872 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
499 :
500 872 : CALL timeset(routineN, handle)
501 :
502 872 : stol = 1.0E-8_dp
503 872 : ssize = 0.2_dp
504 872 : maxit = 999
505 872 : fail = .FALSE.
506 872 : bisec = .FALSE.
507 :
508 105260 : dg = 0._dp
509 :
510 : CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
511 872 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
512 :
513 61676 : DO i = 1, nrow_local
514 60804 : j = row_indices(i)
515 16256708 : DO k = 1, ncol_local
516 16195032 : l = col_indices(k)
517 16255836 : dg(l) = dg(l) + local_data(i, k)*g(j)
518 : END DO
519 : END DO
520 872 : CALL para_env%sum(dg)
521 :
522 872 : set = .FALSE.
523 :
524 : DO
525 :
526 : ! calculating Lambda
527 :
528 872 : lp = 0.0_dp
529 872 : iref = 1
530 872 : ln = 0.0_dp
531 872 : IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
532 :
533 872 : iter = 0
534 : DO
535 2565 : iter = iter + 1
536 2565 : fun = 0.0_dp
537 2565 : fung = 0.0_dp
538 264159 : DO indf = 1, ndf
539 261594 : fun = fun + dg(indf)**2/(ln - eigval(indf))
540 264159 : fung = fung - dg(indf)**2/(ln - eigval(indf)**2)
541 : END DO
542 2565 : fun = fun - ln
543 2565 : fung = fung - one
544 2565 : step = fun/fung
545 2565 : ln = ln - step
546 2565 : IF (ABS(step) < stol) GOTO 200
547 1694 : IF (iter >= maxit) EXIT
548 : END DO
549 : 100 CONTINUE
550 1 : bisec = .TRUE.
551 1 : iter = 0
552 1 : maxit = 9999
553 1 : lam1 = 0.0_dp
554 1 : IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
555 : fun1 = 0.0_dp
556 31 : DO indf = 1, ndf
557 31 : fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
558 : END DO
559 1 : fun1 = fun1 - lam1
560 1 : step = ABS(lam1)/1000.0_dp
561 1 : IF (step < ssize) step = ssize
562 : DO
563 1 : iter = iter + 1
564 1 : IF (iter > maxit) THEN
565 : ln = 0.0_dp
566 872 : lp = 0.0_dp
567 : fail = .TRUE.
568 : GOTO 300
569 : END IF
570 1 : fun2 = 0.0_dp
571 1 : lam2 = lam1 - iter*step
572 31 : DO indf = 1, ndf
573 31 : fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf))
574 : END DO
575 1 : fun2 = fun2 - lam2
576 873 : IF (fun2*fun1 < 0.0_dp) THEN
577 : iter = 0
578 : DO
579 25 : iter = iter + 1
580 25 : IF (iter > maxit) THEN
581 : ln = 0.0_dp
582 : lp = 0.0_dp
583 : fail = .TRUE.
584 : GOTO 300
585 : END IF
586 25 : step = (lam1 + lam2)/2
587 25 : fun3 = 0.0_dp
588 775 : DO indf = 1, ndf
589 775 : fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
590 : END DO
591 25 : fun3 = fun3 - step
592 :
593 25 : IF (ABS(step - lam2) < stol) THEN
594 : ln = step
595 : GOTO 200
596 : END IF
597 :
598 24 : IF (fun3*fun1 < stol) THEN
599 : lam2 = step
600 : ELSE
601 24 : lam1 = step
602 : END IF
603 : END DO
604 : END IF
605 : END DO
606 :
607 : 200 CONTINUE
608 873 : IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
609 : (eigval(iref) > 0.0_dp))) THEN
610 :
611 1 : IF (.NOT. bisec) GOTO 100
612 : ln = 0.0_dp
613 : lp = 0.0_dp
614 : fail = .TRUE.
615 : END IF
616 :
617 : 300 CONTINUE
618 :
619 872 : IF (fail .AND. .NOT. set) THEN
620 0 : set = .TRUE.
621 0 : DO indf = 1, ndf
622 0 : eigval(indf) = eigval(indf)*work(indf)
623 : END DO
624 : CYCLE
625 : END IF
626 :
627 872 : IF (.NOT. set) THEN
628 105260 : work(1:ndf) = one
629 : END IF
630 :
631 105260 : DO indf = 1, ndf
632 105260 : eigval(indf) = eigval(indf) - ln
633 : END DO
634 : EXIT
635 : END DO
636 :
637 872 : CALL timestop(handle)
638 :
639 872 : END SUBROUTINE rat_fun_opt
640 :
641 : ! **************************************************************************************************
642 : !> \brief ...
643 : !> \param ndf ...
644 : !> \param dx ...
645 : !> \param dg ...
646 : !> \param hess_mat ...
647 : !> \param work ...
648 : !> \param para_env ...
649 : ! **************************************************************************************************
650 10640 : SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env)
651 : INTEGER, INTENT(IN) :: ndf
652 : REAL(KIND=dp), INTENT(INOUT) :: dx(ndf), dg(ndf)
653 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
654 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
655 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
656 :
657 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bfgs'
658 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
659 :
660 : INTEGER :: handle, i, j, k, l, ncol_local, &
661 : nrow_local
662 5320 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
663 : REAL(KIND=dp) :: DDOT, dxw, gdx
664 5320 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_hes
665 :
666 5320 : CALL timeset(routineN, handle)
667 :
668 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
669 5320 : local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
670 :
671 582085 : work = zero
672 316906 : DO i = 1, nrow_local
673 311586 : j = row_indices(i)
674 53992366 : DO k = 1, ncol_local
675 53675460 : l = col_indices(k)
676 53987046 : work(j) = work(j) + local_hes(i, k)*dx(l)
677 : END DO
678 : END DO
679 :
680 5320 : CALL para_env%sum(work)
681 :
682 5320 : gdx = DDOT(ndf, dg, 1, dx, 1)
683 5320 : gdx = one/gdx
684 5320 : dxw = DDOT(ndf, dx, 1, work, 1)
685 5320 : dxw = one/dxw
686 :
687 316906 : DO i = 1, nrow_local
688 311586 : j = row_indices(i)
689 53992366 : DO k = 1, ncol_local
690 53675460 : l = col_indices(k)
691 : local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
692 53987046 : dxw*work(j)*work(l)
693 : END DO
694 : END DO
695 :
696 5320 : CALL timestop(handle)
697 :
698 5320 : END SUBROUTINE bfgs
699 :
700 : ! **************************************************************************************************
701 : !> \brief ...
702 : !> \param ndf ...
703 : !> \param eigval ...
704 : !> \param work ...
705 : ! **************************************************************************************************
706 872 : SUBROUTINE set_hes_eig(ndf, eigval, work)
707 : INTEGER, INTENT(IN) :: ndf
708 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf), work(ndf)
709 :
710 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_hes_eig'
711 : REAL(KIND=dp), PARAMETER :: max_neg = -0.5_dp, max_pos = 5.0_dp, &
712 : min_eig = 0.005_dp, one = 1.0_dp
713 :
714 : INTEGER :: handle, indf
715 : LOGICAL :: neg
716 :
717 872 : CALL timeset(routineN, handle)
718 :
719 105260 : DO indf = 1, ndf
720 104388 : IF (eigval(indf) < 0.0_dp) neg = .TRUE.
721 105260 : IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
722 : END DO
723 105260 : DO indf = 1, ndf
724 105260 : IF (eigval(indf) < 0.0_dp) THEN
725 1 : IF (eigval(indf) < max_neg) THEN
726 0 : eigval(indf) = max_neg
727 1 : ELSE IF (eigval(indf) > -min_eig) THEN
728 1 : eigval(indf) = -min_eig
729 : END IF
730 104387 : ELSE IF (eigval(indf) < 1000.0_dp) THEN
731 104387 : IF (eigval(indf) < min_eig) THEN
732 272 : eigval(indf) = min_eig
733 104115 : ELSE IF (eigval(indf) > max_pos) THEN
734 0 : eigval(indf) = max_pos
735 : END IF
736 : END IF
737 : END DO
738 :
739 105260 : DO indf = 1, ndf
740 105260 : IF (eigval(indf) < 0.0_dp) THEN
741 1 : work(indf) = -one
742 : ELSE
743 104387 : work(indf) = one
744 : END IF
745 : END DO
746 :
747 872 : CALL timestop(handle)
748 :
749 872 : END SUBROUTINE set_hes_eig
750 :
751 : ! **************************************************************************************************
752 : !> \brief ...
753 : !> \param ndf ...
754 : !> \param eigval ...
755 : !> \param eigvec_mat ...
756 : !> \param hess_tmp ...
757 : !> \param dr ...
758 : !> \param g ...
759 : !> \param para_env ...
760 : !> \param use_rfo ...
761 : ! **************************************************************************************************
762 19647 : SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
763 :
764 : INTEGER, INTENT(IN) :: ndf
765 : REAL(KIND=dp), INTENT(INOUT) :: eigval(ndf)
766 : TYPE(cp_fm_type), INTENT(IN) :: eigvec_mat, hess_tmp
767 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
768 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
769 : LOGICAL :: use_rfo
770 :
771 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
772 :
773 : INTEGER :: i, indf, j, k, l, ncol_local, nrow_local
774 6549 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
775 6549 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
776 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
777 : TYPE(cp_fm_type) :: tmp
778 :
779 6549 : CALL cp_fm_to_fm(eigvec_mat, hess_tmp)
780 6549 : IF (use_rfo) THEN
781 105260 : DO indf = 1, ndf
782 105260 : eigval(indf) = one/eigval(indf)
783 : END DO
784 : ELSE
785 565450 : DO indf = 1, ndf
786 565450 : eigval(indf) = one/MAX(0.0001_dp, eigval(indf))
787 : END DO
788 : END IF
789 :
790 6549 : CALL cp_fm_column_scale(hess_tmp, eigval)
791 6549 : CALL cp_fm_get_info(eigvec_mat, matrix_struct=matrix_struct)
792 6549 : CALL cp_fm_create(tmp, matrix_struct, name="tmp")
793 :
794 6549 : CALL parallel_gemm("N", "T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
795 :
796 6549 : CALL cp_fm_transpose(tmp, hess_tmp)
797 6549 : CALL cp_fm_release(tmp)
798 :
799 : ! ** New step **
800 :
801 : CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
802 6549 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
803 :
804 670710 : dr = 0.0_dp
805 363738 : DO i = 1, nrow_local
806 357189 : j = row_indices(i)
807 63787395 : DO k = 1, ncol_local
808 63423657 : l = col_indices(k)
809 63780846 : dr(j) = dr(j) - local_data(i, k)*g(l)
810 : END DO
811 : END DO
812 :
813 6549 : CALL para_env%sum(dr)
814 :
815 6549 : END SUBROUTINE geoopt_get_step
816 :
817 : ! **************************************************************************************************
818 : !> \brief ...
819 : !> \param ndf ...
820 : !> \param step ...
821 : !> \param rad ...
822 : !> \param rat ...
823 : !> \param dr ...
824 : !> \param output_unit ...
825 : ! **************************************************************************************************
826 6549 : SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
827 : INTEGER, INTENT(IN) :: ndf
828 : REAL(KIND=dp), INTENT(INOUT) :: step, rad, rat, dr(ndf)
829 : INTEGER, INTENT(IN) :: output_unit
830 :
831 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trust_radius'
832 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
833 :
834 : INTEGER :: handle
835 : REAL(KIND=dp) :: scal
836 :
837 6549 : CALL timeset(routineN, handle)
838 :
839 677259 : step = MAXVAL(ABS(dr))
840 6549 : scal = MAX(one, rad/step)
841 :
842 6549 : IF (step > rad) THEN
843 436 : rat = rad/step
844 436 : CALL DSCAL(ndf, rat, dr, 1)
845 436 : step = rad
846 436 : IF (output_unit > 0) THEN
847 : WRITE (unit=output_unit, FMT="(/,T2,A,F8.5)") &
848 221 : " Step is scaled; Scaling factor = ", rat
849 221 : CALL m_flush(output_unit)
850 : END IF
851 : END IF
852 6549 : CALL timestop(handle)
853 :
854 6549 : END SUBROUTINE trust_radius
855 :
856 : ! **************************************************************************************************
857 : !> \brief ...
858 : !> \param ndf ...
859 : !> \param work ...
860 : !> \param hess_mat ...
861 : !> \param dr ...
862 : !> \param g ...
863 : !> \param conv ...
864 : !> \param pred ...
865 : !> \param para_env ...
866 : ! **************************************************************************************************
867 13098 : SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
868 :
869 : INTEGER, INTENT(IN) :: ndf
870 : REAL(KIND=dp), INTENT(INOUT) :: work(ndf)
871 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
872 : REAL(KIND=dp), INTENT(INOUT) :: dr(ndf), g(ndf)
873 : LOGICAL, INTENT(INOUT) :: conv
874 : REAL(KIND=dp), INTENT(INOUT) :: pred
875 : TYPE(mp_para_env_type), POINTER :: para_env
876 :
877 : CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_predict'
878 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
879 :
880 : INTEGER :: handle, i, j, k, l, ncol_local, &
881 : nrow_local
882 6549 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
883 : REAL(KIND=dp) :: DDOT, ener1, ener2
884 6549 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: local_data
885 :
886 6549 : CALL timeset(routineN, handle)
887 :
888 6549 : ener1 = DDOT(ndf, g, 1, dr, 1)
889 :
890 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
891 6549 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
892 :
893 670710 : work = zero
894 363738 : DO i = 1, nrow_local
895 357189 : j = row_indices(i)
896 63787395 : DO k = 1, ncol_local
897 63423657 : l = col_indices(k)
898 63780846 : work(j) = work(j) + local_data(i, k)*dr(l)
899 : END DO
900 : END DO
901 :
902 6549 : CALL para_env%sum(work)
903 6549 : ener2 = DDOT(ndf, dr, 1, work, 1)
904 6549 : pred = ener1 + 0.5_dp*ener2
905 6549 : conv = .FALSE.
906 6549 : CALL timestop(handle)
907 :
908 6549 : END SUBROUTINE energy_predict
909 :
910 : ! **************************************************************************************************
911 : !> \brief ...
912 : !> \param rat ...
913 : !> \param rad ...
914 : !> \param step ...
915 : !> \param ediff ...
916 : ! **************************************************************************************************
917 777 : SUBROUTINE update_trust_rad(rat, rad, step, ediff)
918 :
919 : REAL(KIND=dp), INTENT(INOUT) :: rat, rad, step, ediff
920 :
921 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_trust_rad'
922 : REAL(KIND=dp), PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp
923 :
924 : INTEGER :: handle
925 :
926 777 : CALL timeset(routineN, handle)
927 :
928 777 : IF (rat > 4.0_dp) THEN
929 0 : IF (ediff < 0.0_dp) THEN
930 0 : rad = step*0.5_dp
931 : ELSE
932 0 : rad = step*0.25_dp
933 : END IF
934 777 : ELSE IF (rat > 2.0_dp) THEN
935 0 : IF (ediff < 0.0_dp) THEN
936 0 : rad = step*0.75_dp
937 : ELSE
938 0 : rad = step*0.5_dp
939 : END IF
940 777 : ELSE IF (rat > 4.0_dp/3.0_dp) THEN
941 0 : IF (ediff < 0.0_dp) THEN
942 0 : rad = step
943 : ELSE
944 0 : rad = step*0.75_dp
945 : END IF
946 777 : ELSE IF (rat > 10.0_dp/9.0_dp) THEN
947 0 : IF (ediff < 0.0_dp) THEN
948 0 : rad = step*1.25_dp
949 : ELSE
950 0 : rad = step
951 : END IF
952 777 : ELSE IF (rat > 0.9_dp) THEN
953 60 : IF (ediff < 0.0_dp) THEN
954 60 : rad = step*1.5_dp
955 : ELSE
956 0 : rad = step*1.25_dp
957 : END IF
958 717 : ELSE IF (rat > 0.75_dp) THEN
959 116 : IF (ediff < 0.0_dp) THEN
960 113 : rad = step*1.25_dp
961 : ELSE
962 3 : rad = step
963 : END IF
964 601 : ELSE IF (rat > 0.5_dp) THEN
965 93 : IF (ediff < 0.0_dp) THEN
966 93 : rad = step
967 : ELSE
968 0 : rad = step*0.75_dp
969 : END IF
970 508 : ELSE IF (rat > 0.25_dp) THEN
971 5 : IF (ediff < 0.0_dp) THEN
972 5 : rad = step*0.75_dp
973 : ELSE
974 0 : rad = step*0.5_dp
975 : END IF
976 503 : ELSE IF (ediff < 0.0_dp) THEN
977 503 : rad = step*0.5_dp
978 : ELSE
979 0 : rad = step*0.25_dp
980 : END IF
981 :
982 777 : rad = MAX(rad, min_trust)
983 777 : rad = MIN(rad, max_trust)
984 777 : CALL timestop(handle)
985 :
986 777 : END SUBROUTINE update_trust_rad
987 :
988 : ! **************************************************************************************************
989 :
990 : ! **************************************************************************************************
991 : !> \brief ...
992 : !> \param geo_section ...
993 : !> \param hess_mat ...
994 : !> \param logger ...
995 : ! **************************************************************************************************
996 5020 : SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger)
997 :
998 : TYPE(section_vals_type), POINTER :: geo_section
999 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1000 : TYPE(cp_logger_type), POINTER :: logger
1001 :
1002 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bfgs_hessian'
1003 :
1004 : INTEGER :: handle, hesunit
1005 :
1006 5020 : CALL timeset(routineN, handle)
1007 :
1008 : hesunit = cp_print_key_unit_nr(logger, geo_section, "BFGS%RESTART", &
1009 : extension=".Hessian", file_form="UNFORMATTED", file_action="WRITE", &
1010 5020 : file_position="REWIND")
1011 :
1012 5020 : CALL cp_fm_write_unformatted(hess_mat, hesunit)
1013 :
1014 5020 : CALL cp_print_key_finished_output(hesunit, logger, geo_section, "BFGS%RESTART")
1015 :
1016 5020 : CALL timestop(handle)
1017 :
1018 5020 : END SUBROUTINE write_bfgs_hessian
1019 : ! **************************************************************************************************
1020 :
1021 : ! **************************************************************************************************
1022 : !> \brief ...
1023 : !> \param force_env ...
1024 : !> \param hess_mat ...
1025 : ! **************************************************************************************************
1026 1023 : SUBROUTINE construct_initial_hess(force_env, hess_mat)
1027 :
1028 : TYPE(force_env_type), POINTER :: force_env
1029 : TYPE(cp_fm_type), INTENT(IN) :: hess_mat
1030 :
1031 : INTEGER :: i, iat_col, iat_row, iglobal, iind, j, &
1032 : jat_row, jglobal, jind, k, natom, &
1033 : ncol_local, nrow_local, z
1034 1023 : INTEGER, ALLOCATABLE, DIMENSION(:) :: at_row
1035 1023 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1036 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_ij, rho_ij
1037 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_ij
1038 : REAL(KIND=dp), DIMENSION(3, 3) :: alpha, r0
1039 1023 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), POINTER :: fixed, local_data
1040 : TYPE(cell_type), POINTER :: cell
1041 : TYPE(cp_subsys_type), POINTER :: subsys
1042 : TYPE(particle_list_type), POINTER :: particles
1043 :
1044 1023 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1045 : CALL cp_subsys_get(subsys, &
1046 1023 : particles=particles)
1047 :
1048 4092 : alpha(1, :) = [1._dp, 0.3949_dp, 0.3949_dp]
1049 4092 : alpha(2, :) = [0.3494_dp, 0.2800_dp, 0.2800_dp]
1050 4092 : alpha(3, :) = [0.3494_dp, 0.2800_dp, 0.1800_dp]
1051 :
1052 4092 : r0(1, :) = [1.35_dp, 2.10_dp, 2.53_dp]
1053 4092 : r0(2, :) = [2.10_dp, 2.87_dp, 3.40_dp]
1054 4092 : r0(3, :) = [2.53_dp, 3.40_dp, 3.40_dp]
1055 :
1056 : CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1057 1023 : local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1058 1023 : natom = particles%n_els
1059 3069 : ALLOCATE (at_row(natom))
1060 4092 : ALLOCATE (rho_ij(natom, natom))
1061 3069 : ALLOCATE (d_ij(natom, natom))
1062 5115 : ALLOCATE (r_ij(natom, natom, 3))
1063 3069 : ALLOCATE (fixed(3, natom))
1064 72775 : fixed = 1.0_dp
1065 1023 : CALL fix_atom_control(force_env, fixed)
1066 4092 : DO i = 1, 3
1067 111720 : CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1068 : END DO
1069 916517 : rho_ij = 0
1070 : !XXXX insert proper rows !XXX
1071 18961 : at_row = 3
1072 18961 : DO i = 1, natom
1073 17938 : CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, z=z)
1074 17938 : IF (z <= 10) at_row(i) = 2
1075 36899 : IF (z <= 2) at_row(i) = 1
1076 : END DO
1077 17938 : DO i = 2, natom
1078 16915 : iat_row = at_row(i)
1079 457747 : DO j = 1, i - 1
1080 439809 : jat_row = at_row(j)
1081 : !pbc for a distance vector
1082 1759236 : r_ij(j, i, :) = pbc(particles%els(i)%r, particles%els(j)%r, cell)
1083 1759236 : r_ij(i, j, :) = -r_ij(j, i, :)
1084 1759236 : d_ij(j, i) = SQRT(DOT_PRODUCT(r_ij(j, i, :), r_ij(j, i, :)))
1085 439809 : d_ij(i, j) = d_ij(j, i)
1086 439809 : rho_ij(j, i) = EXP(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1087 456724 : rho_ij(i, j) = rho_ij(j, i)
1088 : END DO
1089 : END DO
1090 54837 : DO i = 1, ncol_local
1091 53814 : iglobal = col_indices(i)
1092 53814 : iind = MOD(iglobal - 1, 3) + 1
1093 53814 : iat_col = (iglobal + 2)/3
1094 53814 : IF (iat_col > natom) CYCLE
1095 4183119 : DO j = 1, nrow_local
1096 4128282 : jglobal = row_indices(j)
1097 4128282 : jind = MOD(jglobal - 1, 3) + 1
1098 4128282 : iat_row = (jglobal + 2)/3
1099 4128282 : IF (iat_row > natom) CYCLE
1100 4128282 : IF (iat_row /= iat_col) THEN
1101 4041846 : IF (d_ij(iat_row, iat_col) < 6.0_dp) &
1102 : local_data(j, i) = local_data(j, i) + &
1103 344520 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1104 : ELSE
1105 : local_data(j, i) = local_data(j, i) + &
1106 86436 : angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1107 : END IF
1108 4128282 : IF (iat_col /= iat_row) THEN
1109 4041846 : IF (d_ij(iat_row, iat_col) < 6.0_dp) &
1110 : local_data(j, i) = local_data(j, i) - &
1111 : dist_second_deriv(r_ij(iat_col, iat_row, :), &
1112 2411640 : iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1113 : ELSE
1114 4214718 : DO k = 1, natom
1115 4128282 : IF (k == iat_col) CYCLE
1116 4041846 : IF (d_ij(iat_row, k) < 6.0_dp) &
1117 : local_data(j, i) = local_data(j, i) + &
1118 : dist_second_deriv(r_ij(iat_col, k, :), &
1119 2498076 : iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1120 : END DO
1121 : END IF
1122 4182096 : IF (fixed(jind, iat_row) < 0.5_dp .OR. fixed(iind, iat_col) < 0.5_dp) THEN
1123 10161 : local_data(j, i) = 0.0_dp
1124 10161 : IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
1125 : END IF
1126 : END DO
1127 : END DO
1128 1023 : DEALLOCATE (fixed)
1129 1023 : DEALLOCATE (rho_ij)
1130 1023 : DEALLOCATE (d_ij)
1131 1023 : DEALLOCATE (r_ij)
1132 1023 : DEALLOCATE (at_row)
1133 :
1134 2046 : END SUBROUTINE construct_initial_hess
1135 :
1136 : ! **************************************************************************************************
1137 : !> \brief ...
1138 : !> \param r1 ...
1139 : !> \param i ...
1140 : !> \param j ...
1141 : !> \param d ...
1142 : !> \param rho ...
1143 : !> \return ...
1144 : ! **************************************************************************************************
1145 689040 : FUNCTION dist_second_deriv(r1, i, j, d, rho) RESULT(deriv)
1146 : REAL(KIND=dp), DIMENSION(3) :: r1
1147 : INTEGER :: i, j
1148 : REAL(KIND=dp) :: d, rho, deriv
1149 :
1150 689040 : deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
1151 689040 : END FUNCTION dist_second_deriv
1152 :
1153 : ! **************************************************************************************************
1154 : !> \brief ...
1155 : !> \param r_ij ...
1156 : !> \param d_ij ...
1157 : !> \param rho_ij ...
1158 : !> \param idir ...
1159 : !> \param jdir ...
1160 : !> \param iat_der ...
1161 : !> \param jat_der ...
1162 : !> \param natom ...
1163 : !> \return ...
1164 : ! **************************************************************************************************
1165 430956 : FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom) RESULT(deriv)
1166 : REAL(KIND=dp), DIMENSION(:, :, :) :: r_ij
1167 : REAL(KIND=dp), DIMENSION(:, :) :: d_ij, rho_ij
1168 : INTEGER :: idir, jdir, iat_der, jat_der, natom
1169 : REAL(KIND=dp) :: deriv
1170 :
1171 : INTEGER :: i, iat, idr, j, jat, jdr
1172 : REAL(KIND=dp) :: d12, d23, d31, D_mat(3, 2), denom1, &
1173 : denom2, denom3, ka1, ka2, ka3, rho12, &
1174 : rho23, rho31, rsst1, rsst2, rsst3
1175 : REAL(KIND=dp), DIMENSION(3) :: r12, r23, r31
1176 :
1177 430956 : deriv = 0._dp
1178 430956 : IF (iat_der == jat_der) THEN
1179 4128282 : DO i = 1, natom - 1
1180 4041846 : IF (rho_ij(iat_der, i) < 0.00001) CYCLE
1181 25705161 : DO j = i + 1, natom
1182 24678540 : IF (rho_ij(iat_der, j) < 0.00001) CYCLE
1183 6052608 : IF (i == iat_der .OR. j == iat_der) CYCLE
1184 6052608 : IF (iat_der < i .OR. iat_der > j) THEN
1185 38921670 : r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1186 3892167 : d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1187 3892167 : rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1188 : ELSE
1189 21604410 : r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1190 2160441 : d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1191 2160441 : rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1192 : END IF
1193 6052608 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1194 60526080 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1195 6052608 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1196 6052608 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1197 6052608 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1198 6052608 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1199 6052608 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1200 6052608 : D_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1201 6052608 : D_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1202 6052608 : D_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1203 6052608 : D_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
1204 : D_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
1205 6052608 : rsst3*r12(idir)/(d31*d12**3)
1206 : D_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1207 6052608 : rsst3*r12(jdir)/(d31*d12**3)
1208 6052608 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1209 6052608 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1210 6052608 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1211 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1212 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1213 28720386 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1214 :
1215 : END DO
1216 : END DO
1217 : ELSE
1218 10312416 : DO i = 1, natom
1219 9967896 : IF (i == iat_der .OR. i == jat_der) CYCLE
1220 9278856 : IF (jat_der < iat_der) THEN
1221 4639428 : iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1222 : ELSE
1223 4639428 : iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1224 : END IF
1225 9278856 : IF (jat < i .OR. iat > i) THEN
1226 70110900 : r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1227 7011090 : d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1228 7011090 : rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1229 : ELSE
1230 22677660 : r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1231 2267766 : d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1232 2267766 : rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1233 : END IF
1234 9278856 : ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1235 92788560 : rsst1 = DOT_PRODUCT(r12, r23); rsst2 = DOT_PRODUCT(r23, r31); rsst3 = DOT_PRODUCT(r31, r12)
1236 9278856 : denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1237 9278856 : denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1238 9278856 : denom1 = SIGN(1.0_dp, denom1)*MAX(ABS(denom1), 0.01_dp)
1239 9278856 : denom2 = SIGN(1.0_dp, denom2)*MAX(ABS(denom2), 0.01_dp)
1240 9278856 : denom3 = SIGN(1.0_dp, denom3)*MAX(ABS(denom3), 0.01_dp)
1241 9278856 : D_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1242 9278856 : D_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
1243 : D_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
1244 9278856 : rsst3*r12(idr)/(d31*d12**3)
1245 9278856 : IF (jat < i .OR. iat > i) THEN
1246 : D_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
1247 7011090 : rsst1*r23(jdr)/(d12*d23**3)
1248 7011090 : D_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1249 7011090 : D_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1250 : ELSE
1251 2267766 : D_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
1252 : D_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
1253 2267766 : rsst2*r31(jdr)/(d23*d31**3)
1254 2267766 : D_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1255 : END IF
1256 9278856 : IF (ABS(denom1) <= 0.011_dp) D_mat(1, 1) = 0.0_dp
1257 9278856 : IF (ABS(denom2) <= 0.011_dp) D_mat(2, 1) = 0.0_dp
1258 9278856 : IF (ABS(denom3) <= 0.011_dp) D_mat(3, 1) = 0.0_dp
1259 :
1260 : deriv = deriv + ka1*D_mat(1, 1)*D_mat(1, 2)/denom1 + &
1261 : ka2*D_mat(2, 1)*D_mat(2, 2)/denom2 + &
1262 10312416 : ka3*D_mat(3, 1)*D_mat(3, 2)/denom3
1263 : END DO
1264 : END IF
1265 430956 : deriv = 0.25_dp*deriv
1266 :
1267 430956 : END FUNCTION angle_second_deriv
1268 :
1269 : END MODULE bfgs_optimizer
|
