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