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