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