Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 956 : 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 478 : CALL timeset(routineN, handle)
99 :
100 478 : NULLIFY (spgr)
101 478 : logger => cp_get_default_logger()
102 478 : spgr => gopt_env%spgr
103 :
104 : output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
105 478 : extension=".geoLog")
106 478 : CALL print_geo_opt_header(gopt_env, output_unit, "CONJUGATE GRADIENTS")
107 :
108 : ! find space_group
109 478 : CALL force_env_get(force_env, subsys=subsys)
110 478 : CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
111 478 : 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 478 : gopt_env, do_update=my_do_update)
139 :
140 : CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
141 478 : "PRINT%PROGRAM_RUN_INFO")
142 478 : IF (PRESENT(do_update)) do_update = my_do_update
143 :
144 478 : CALL timestop(handle)
145 :
146 478 : END SUBROUTINE geoopt_cg
147 :
148 : ! **************************************************************************************************
149 : !> \brief This really performs the conjugate gradients optimization
150 : !> \param force_env ...
151 : !> \param x0 ...
152 : !> \param gopt_param ...
153 : !> \param output_unit ...
154 : !> \param globenv ...
155 : !> \param gopt_env ...
156 : !> \param do_update ...
157 : !> \par History
158 : !> 10.2005 created [tlaino]
159 : !> \author Teodoro Laino
160 : ! **************************************************************************************************
161 478 : RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
162 : gopt_env, do_update)
163 : TYPE(force_env_type), POINTER :: force_env
164 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
165 : TYPE(gopt_param_type), POINTER :: gopt_param
166 : INTEGER, INTENT(IN) :: output_unit
167 : TYPE(global_environment_type), POINTER :: globenv
168 : TYPE(gopt_f_type), POINTER :: gopt_env
169 : LOGICAL, INTENT(OUT), OPTIONAL :: do_update
170 :
171 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cg_main'
172 :
173 : CHARACTER(LEN=5) :: wildcard
174 : INTEGER :: handle, iter_nr, its, max_steep_steps, &
175 : maxiter
176 : LOGICAL :: conv, Fletcher_Reeves, &
177 : save_consistent_energy_force, &
178 : should_stop
179 : REAL(KIND=dp) :: emin, eold, opt_energy, res_lim, t_diff, &
180 : t_now, t_old
181 478 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xold
182 478 : REAL(KIND=dp), DIMENSION(:), POINTER :: g, h, xi
183 : TYPE(cell_type), POINTER :: cell
184 : TYPE(cp_logger_type), POINTER :: logger
185 : TYPE(cp_subsys_type), POINTER :: subsys
186 : TYPE(section_vals_type), POINTER :: root_section
187 : TYPE(spgr_type), POINTER :: spgr
188 :
189 478 : CALL timeset(routineN, handle)
190 478 : t_old = m_walltime()
191 478 : NULLIFY (logger, g, h, xi, spgr)
192 478 : root_section => force_env%root_section
193 478 : logger => cp_get_default_logger()
194 478 : conv = .FALSE.
195 478 : maxiter = gopt_param%max_iter
196 478 : max_steep_steps = gopt_param%max_steep_steps
197 478 : Fletcher_Reeves = gopt_param%Fletcher_Reeves
198 478 : res_lim = gopt_param%restart_limit
199 1434 : ALLOCATE (g(SIZE(x0)))
200 1434 : ALLOCATE (h(SIZE(x0)))
201 1434 : ALLOCATE (xi(SIZE(x0)))
202 1434 : ALLOCATE (xold(SIZE(x0)))
203 478 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
204 :
205 478 : spgr => gopt_env%spgr
206 : ! applies rotation matrices to coordinates
207 478 : IF (spgr%keep_space_group) THEN
208 2 : CALL spgr_apply_rotations_coord(spgr, x0)
209 : END IF
210 :
211 : ! Evaluate energy and forces at the first step
212 : ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
213 478 : save_consistent_energy_force = gopt_env%require_consistent_energy_force
214 478 : gopt_env%require_consistent_energy_force = .FALSE.
215 :
216 : CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
217 478 : .FALSE., gopt_env%force_env%para_env)
218 :
219 478 : gopt_env%require_consistent_energy_force = save_consistent_energy_force
220 :
221 : ! Symmetrize coordinates and forces
222 478 : IF (spgr%keep_space_group) THEN
223 2 : CALL spgr_apply_rotations_coord(spgr, x0)
224 2 : CALL spgr_apply_rotations_force(spgr, xi)
225 : END IF
226 :
227 185012 : g = -xi
228 185012 : h = g
229 185012 : xi = h
230 478 : emin = HUGE(0.0_dp)
231 478 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
232 : ! Main Loop
233 478 : wildcard = " SD"
234 478 : t_now = m_walltime()
235 478 : t_diff = t_now - t_old
236 478 : t_old = t_now
237 478 : CALL gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, used_time=t_diff, its=iter_nr)
238 478 : eold = opt_energy
239 2526 : DO its = iter_nr + 1, maxiter
240 2526 : CALL cp_iterate(logger%iter_info, last=(its == maxiter))
241 2526 : CALL section_vals_val_set(gopt_env%geo_section, "STEP_START_VAL", i_val=its)
242 2526 : CALL gopt_f_ii(its, output_unit)
243 :
244 : ! Symmetrize coordinates and forces
245 2526 : IF (spgr%keep_space_group) THEN
246 2 : CALL spgr_apply_rotations_coord(spgr, x0)
247 2 : CALL spgr_apply_rotations_force(spgr, g)
248 2 : CALL spgr_apply_rotations_force(spgr, xi)
249 : END IF
250 :
251 415038 : xold(:) = x0
252 :
253 : ! Line minimization
254 2526 : CALL cg_linmin(gopt_env, x0, xi, g, opt_energy, output_unit, gopt_param, globenv)
255 :
256 : ! Applies rotation matrices to coordinates
257 2526 : IF (spgr%keep_space_group) THEN
258 2 : CALL spgr_apply_rotations_coord(spgr, x0)
259 : END IF
260 :
261 : ! Check for an external exit command
262 2526 : CALL external_control(should_stop, "GEO", globenv=globenv)
263 2526 : IF (should_stop) EXIT
264 :
265 : ! Some IO and Convergence check
266 2526 : t_now = m_walltime()
267 2526 : t_diff = t_now - t_old
268 2526 : t_old = t_now
269 : CALL gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
270 : output_unit, eold, emin, wildcard, gopt_param, SIZE(x0), x0 - xold, xi, conv, &
271 415038 : used_time=t_diff)
272 2526 : eold = opt_energy
273 2526 : emin = MIN(emin, opt_energy)
274 :
275 2526 : IF (conv .OR. (its == maxiter)) EXIT
276 : ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
277 2048 : save_consistent_energy_force = gopt_env%require_consistent_energy_force
278 2048 : gopt_env%require_consistent_energy_force = .FALSE.
279 :
280 : CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
281 2048 : .FALSE., gopt_env%force_env%para_env)
282 :
283 2048 : gopt_env%require_consistent_energy_force = save_consistent_energy_force
284 :
285 : ! Symmetrize coordinates and forces
286 2048 : IF (spgr%keep_space_group) THEN
287 0 : CALL spgr_apply_rotations_force(spgr, xi)
288 : END IF
289 :
290 : ! Get Conjugate Directions: updates the searching direction (h)
291 2048 : wildcard = " CG"
292 2048 : CALL get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
293 :
294 : ! Symmetrize coordinates and forces
295 2048 : IF (spgr%keep_space_group) THEN
296 0 : CALL spgr_apply_rotations_force(spgr, g)
297 0 : CALL spgr_apply_rotations_force(spgr, h)
298 : END IF
299 :
300 : ! Reset Condition or Steepest Descent Requested
301 : ! ABS(DOT_PRODUCT(g, h))/SQRT((DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h))) > res_lim ...
302 : IF ((DOT_PRODUCT(g, h)*DOT_PRODUCT(g, h)) > (res_lim*res_lim*DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h)) &
303 967596 : .OR. its + 1 <= max_steep_steps) THEN
304 : ! Steepest Descent
305 560 : wildcard = " SD"
306 80728 : h = -xi
307 : END IF
308 645064 : g = -xi
309 648068 : xi = h
310 : END DO
311 :
312 478 : IF (its == maxiter .AND. (.NOT. conv)) THEN
313 84 : CALL print_geo_opt_nc(gopt_env, output_unit)
314 : END IF
315 :
316 : ! Write final particle information and restart, if converged
317 478 : IF (PRESENT(do_update)) do_update = conv
318 478 : CALL cp_iterate(logger%iter_info, last=.TRUE., increment=0)
319 : CALL gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
320 478 : gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
321 :
322 478 : DEALLOCATE (xold)
323 478 : DEALLOCATE (g)
324 478 : DEALLOCATE (h)
325 478 : DEALLOCATE (xi)
326 :
327 478 : CALL timestop(handle)
328 :
329 478 : END SUBROUTINE cp_cg_main
330 :
331 : END MODULE cg_optimizer
|
