Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for Geometry optimization using Conjugate Gradients
10 : !> \author Teodoro Laino [teo]
11 : !> 10.2005
12 : ! **************************************************************************************************
13 : MODULE cg_optimizer
14 :
15 : USE cell_types, ONLY: cell_type
16 : USE cg_utils, ONLY: cg_linmin,&
17 : get_conjugate_direction
18 : USE cp_external_control, ONLY: external_control
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_type
21 : USE cp_output_handling, ONLY: cp_iterate,&
22 : cp_print_key_finished_output,&
23 : cp_print_key_unit_nr
24 : USE cp_subsys_types, ONLY: cp_subsys_type
25 : USE force_env_types, ONLY: force_env_get,&
26 : force_env_type
27 : USE global_types, ONLY: global_environment_type
28 : USE gopt_f_methods, ONLY: cp_eval_at,&
29 : gopt_f_ii,&
30 : gopt_f_io,&
31 : gopt_f_io_finalize,&
32 : gopt_f_io_init,&
33 : print_geo_opt_header,&
34 : print_geo_opt_nc
35 : USE gopt_f_types, ONLY: gopt_f_type
36 : USE gopt_param_types, ONLY: gopt_param_type
37 : USE input_constants, ONLY: default_cell_method_id,&
38 : default_minimization_method_id,&
39 : default_ts_method_id
40 : USE input_section_types, ONLY: section_vals_type,&
41 : section_vals_val_get,&
42 : section_vals_val_set
43 : USE kinds, ONLY: dp
44 : USE machine, ONLY: m_walltime
45 : USE space_groups, ONLY: identify_space_group,&
46 : print_spgr,&
47 : spgr_apply_rotations_coord,&
48 : spgr_apply_rotations_force
49 : USE space_groups_types, ONLY: spgr_type
50 : #include "../base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 : PRIVATE
54 :
55 : PUBLIC :: geoopt_cg
56 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_optimizer'
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Driver for conjugate gradient optimization technique
63 : !> \param force_env ...
64 : !> \param gopt_param ...
65 : !> \param globenv ...
66 : !> \param geo_section ...
67 : !> \param gopt_env ...
68 : !> \param x0 ...
69 : !> \param do_update ...
70 : !> \par History
71 : !> 10.2005 created [tlaino]
72 : !> \author Teodoro Laino
73 : ! **************************************************************************************************
74 516 : RECURSIVE SUBROUTINE geoopt_cg(force_env, gopt_param, globenv, geo_section, &
75 : gopt_env, x0, do_update)
76 :
77 : TYPE(force_env_type), POINTER :: force_env
78 : TYPE(gopt_param_type), POINTER :: gopt_param
79 : TYPE(global_environment_type), POINTER :: globenv
80 : TYPE(section_vals_type), POINTER :: geo_section
81 : TYPE(gopt_f_type), POINTER :: gopt_env
82 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
83 : LOGICAL, INTENT(OUT), OPTIONAL :: do_update
84 :
85 : CHARACTER(len=*), PARAMETER :: routineN = 'geoopt_cg'
86 :
87 : INTEGER :: handle, output_unit
88 : LOGICAL :: my_do_update
89 : TYPE(cp_logger_type), POINTER :: logger
90 : TYPE(cp_subsys_type), POINTER :: subsys
91 : TYPE(spgr_type), POINTER :: spgr
92 :
93 258 : CALL timeset(routineN, handle)
94 :
95 258 : NULLIFY (spgr)
96 258 : logger => cp_get_default_logger()
97 258 : spgr => gopt_env%spgr
98 :
99 : output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
100 258 : extension=".geoLog")
101 258 : CALL print_geo_opt_header(gopt_env, output_unit, "CONJUGATE GRADIENTS")
102 :
103 : ! find space_group
104 258 : CALL force_env_get(force_env, subsys=subsys)
105 258 : CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
106 258 : IF (spgr%keep_space_group) THEN
107 8 : SELECT CASE (gopt_env%type_id)
108 : CASE (default_minimization_method_id, default_ts_method_id, default_cell_method_id)
109 4 : CALL force_env_get(force_env, subsys=subsys)
110 4 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
111 4 : CALL spgr_apply_rotations_coord(spgr, x0)
112 4 : CALL print_spgr(spgr)
113 : CASE DEFAULT
114 4 : spgr%keep_space_group = .FALSE.
115 : END SELECT
116 : END IF
117 :
118 : CALL cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
119 258 : gopt_env, do_update=my_do_update)
120 :
121 : ! show space_group
122 258 : CALL section_vals_val_get(geo_section, "SHOW_SPACE_GROUP", l_val=spgr%show_space_group)
123 258 : IF (spgr%show_space_group) THEN
124 2 : IF (spgr%keep_space_group) THEN
125 0 : CALL force_env_get(force_env, subsys=subsys)
126 : END IF
127 2 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
128 2 : CALL print_spgr(spgr)
129 : END IF
130 :
131 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
132 258 : "PRINT%PROGRAM_RUN_INFO")
133 258 : IF (PRESENT(do_update)) do_update = my_do_update
134 :
135 258 : CALL timestop(handle)
136 :
137 258 : END SUBROUTINE geoopt_cg
138 :
139 : ! **************************************************************************************************
140 : !> \brief This really performs the conjugate gradients optimization
141 : !> \param force_env ...
142 : !> \param x0 ...
143 : !> \param gopt_param ...
144 : !> \param output_unit ...
145 : !> \param globenv ...
146 : !> \param gopt_env ...
147 : !> \param do_update ...
148 : !> \par History
149 : !> 10.2005 created [tlaino]
150 : !> \author Teodoro Laino
151 : ! **************************************************************************************************
152 258 : RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
153 : gopt_env, do_update)
154 : TYPE(force_env_type), POINTER :: force_env
155 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
156 : TYPE(gopt_param_type), POINTER :: gopt_param
157 : INTEGER, INTENT(IN) :: output_unit
158 : TYPE(global_environment_type), POINTER :: globenv
159 : TYPE(gopt_f_type), POINTER :: gopt_env
160 : LOGICAL, INTENT(OUT), OPTIONAL :: do_update
161 :
162 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cg_main'
163 :
164 : CHARACTER(LEN=5) :: wildcard
165 : INTEGER :: handle, iter_nr, its, max_steep_steps, &
166 : maxiter
167 : LOGICAL :: conv, evaluate_before_io, &
168 : Fletcher_Reeves, &
169 : save_consistent_energy_force, &
170 : should_stop
171 : REAL(KIND=dp) :: emin, eold, opt_energy, res_lim, t_diff, &
172 : t_now, t_old
173 258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xold
174 258 : REAL(KIND=dp), DIMENSION(:), POINTER :: g, h, xi
175 : TYPE(cell_type), POINTER :: cell
176 : TYPE(cp_logger_type), POINTER :: logger
177 : TYPE(cp_subsys_type), POINTER :: subsys
178 : TYPE(section_vals_type), POINTER :: root_section
179 : TYPE(spgr_type), POINTER :: spgr
180 :
181 258 : CALL timeset(routineN, handle)
182 258 : t_old = m_walltime()
183 258 : NULLIFY (logger, g, h, xi, spgr)
184 258 : root_section => force_env%root_section
185 258 : logger => cp_get_default_logger()
186 258 : conv = .FALSE.
187 258 : maxiter = gopt_param%max_iter
188 258 : max_steep_steps = gopt_param%max_steep_steps
189 258 : Fletcher_Reeves = gopt_param%Fletcher_Reeves
190 : ! Dimer evaluations update the rotational state, so retain their original timing.
191 258 : evaluate_before_io = gopt_env%type_id /= default_ts_method_id
192 258 : res_lim = gopt_param%restart_limit
193 774 : ALLOCATE (g(SIZE(x0)))
194 516 : ALLOCATE (h(SIZE(x0)))
195 516 : ALLOCATE (xi(SIZE(x0)))
196 516 : ALLOCATE (xold(SIZE(x0)))
197 258 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
198 :
199 258 : spgr => gopt_env%spgr
200 : ! applies rotation matrices to coordinates
201 258 : IF (spgr%keep_space_group) THEN
202 4 : CALL spgr_apply_rotations_coord(spgr, x0)
203 : END IF
204 :
205 : ! Evaluate energy and forces at the first step
206 : ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
207 258 : save_consistent_energy_force = gopt_env%require_consistent_energy_force
208 258 : gopt_env%require_consistent_energy_force = .FALSE.
209 :
210 : CALL cp_eval_at(gopt_env, x0, opt_energy, xi, master=gopt_env%force_env%para_env%mepos, &
211 258 : para_env=gopt_env%force_env%para_env)
212 :
213 258 : gopt_env%require_consistent_energy_force = save_consistent_energy_force
214 :
215 : ! Symmetrize coordinates and forces
216 258 : IF (spgr%keep_space_group) THEN
217 4 : CALL spgr_apply_rotations_coord(spgr, x0)
218 4 : CALL spgr_apply_rotations_force(spgr, xi)
219 : END IF
220 :
221 144432 : g = -xi
222 144432 : h = g
223 144432 : xi = h
224 258 : emin = HUGE(0.0_dp)
225 258 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
226 : ! Main Loop
227 258 : wildcard = " SD"
228 258 : t_now = m_walltime()
229 258 : t_diff = t_now - t_old
230 258 : t_old = t_now
231 258 : CALL gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, used_time=t_diff, its=iter_nr)
232 258 : eold = opt_energy
233 1848 : DO its = iter_nr + 1, maxiter
234 1848 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
235 1848 : CALL section_vals_val_set(gopt_env%geo_section, "STEP_START_VAL", i_val=its)
236 1848 : CALL gopt_f_ii(its, output_unit)
237 :
238 : ! Symmetrize coordinates and forces
239 1848 : IF (spgr%keep_space_group) THEN
240 66 : CALL spgr_apply_rotations_coord(spgr, x0)
241 66 : CALL spgr_apply_rotations_force(spgr, g)
242 66 : CALL spgr_apply_rotations_force(spgr, xi)
243 : END IF
244 :
245 371178 : xold(:) = x0
246 :
247 : ! Line minimization
248 1848 : CALL cg_linmin(gopt_env, x0, xi, g, opt_energy, output_unit, gopt_param, globenv)
249 :
250 : ! Applies rotation matrices to coordinates
251 1848 : IF (spgr%keep_space_group) THEN
252 66 : CALL spgr_apply_rotations_coord(spgr, x0)
253 : END IF
254 :
255 1848 : IF (evaluate_before_io) THEN
256 : ! Keep the accepted point and the physical state consistent even on external stop.
257 1004 : CALL cg_eval_current(gopt_env, x0, opt_energy, xi, spgr)
258 : END IF
259 :
260 : ! Check for an external exit command
261 1848 : CALL external_control(should_stop, "GEO", globenv=globenv)
262 1848 : IF (should_stop) EXIT
263 :
264 : ! Some IO and Convergence check
265 1848 : t_now = m_walltime()
266 1848 : t_diff = t_now - t_old
267 1848 : t_old = t_now
268 : CALL gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
269 : output_unit, eold, emin, wildcard, gopt_param, SIZE(x0), x0 - xold, xi, conv, &
270 371178 : used_time=t_diff)
271 1848 : eold = opt_energy
272 1848 : emin = MIN(emin, opt_energy)
273 :
274 1848 : IF (conv .OR. (its == maxiter)) EXIT
275 1590 : IF (.NOT. evaluate_before_io) THEN
276 698 : CALL cg_eval_current(gopt_env, x0, opt_energy, xi, spgr)
277 : END IF
278 :
279 : ! Get Conjugate Directions: updates the searching direction (h)
280 1590 : wildcard = " CG"
281 1590 : CALL get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
282 :
283 : ! Symmetrize coordinates and forces
284 1590 : IF (spgr%keep_space_group) THEN
285 62 : CALL spgr_apply_rotations_force(spgr, g)
286 62 : CALL spgr_apply_rotations_force(spgr, h)
287 : END IF
288 :
289 : ! Reset Condition or Steepest Descent Requested
290 : ! ABS(DOT_PRODUCT(g, h))/SQRT((DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h))) > res_lim ...
291 : IF ((DOT_PRODUCT(g, h)*DOT_PRODUCT(g, h)) > (res_lim*res_lim*DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h)) &
292 895296 : .OR. its + 1 <= max_steep_steps) THEN
293 : ! Steepest Descent
294 500 : wildcard = " SD"
295 99424 : h = -xi
296 : END IF
297 597924 : g = -xi
298 600030 : xi = h
299 : END DO
300 :
301 258 : IF (its == maxiter .AND. (.NOT. conv)) THEN
302 76 : CALL print_geo_opt_nc(gopt_env, output_unit)
303 : END IF
304 :
305 : ! Write final particle information and restart, if converged
306 258 : IF (PRESENT(do_update)) do_update = conv
307 258 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
308 : CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
309 258 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
310 :
311 258 : DEALLOCATE (xold)
312 258 : DEALLOCATE (g)
313 258 : DEALLOCATE (h)
314 258 : DEALLOCATE (xi)
315 :
316 258 : CALL timestop(handle)
317 :
318 516 : END SUBROUTINE cp_cg_main
319 :
320 : ! **************************************************************************************************
321 : !> \brief Evaluate the current CG point with the force-consistency setting used by CG.
322 : !> \param gopt_env Optimization environment
323 : !> \param x Current coordinates, possibly projected by cp_eval_at
324 : !> \param energy Energy at the evaluated coordinates
325 : !> \param gradient Gradient at the evaluated coordinates
326 : !> \param spgr Space-group data
327 : ! **************************************************************************************************
328 1702 : SUBROUTINE cg_eval_current(gopt_env, x, energy, gradient, spgr)
329 : TYPE(gopt_f_type), POINTER :: gopt_env
330 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
331 : REAL(KIND=dp), INTENT(OUT) :: energy
332 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
333 : TYPE(spgr_type), POINTER :: spgr
334 :
335 : LOGICAL :: save_consistent_energy_force
336 :
337 1702 : save_consistent_energy_force = gopt_env%require_consistent_energy_force
338 1702 : gopt_env%require_consistent_energy_force = .FALSE.
339 : CALL cp_eval_at(gopt_env, x, energy, gradient, master=gopt_env%force_env%para_env%mepos, &
340 1702 : para_env=gopt_env%force_env%para_env)
341 1702 : gopt_env%require_consistent_energy_force = save_consistent_energy_force
342 :
343 1702 : IF (spgr%keep_space_group) CALL spgr_apply_rotations_force(spgr, gradient)
344 1702 : END SUBROUTINE cg_eval_current
345 :
346 : END MODULE cg_optimizer
|
