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