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 Utilities for Geometry optimization using Conjugate Gradients
10 : !> \author Teodoro Laino [teo]
11 : !> 10.2005
12 : ! **************************************************************************************************
13 : MODULE cg_utils
14 : USE cp_external_control, ONLY: external_control
15 : USE dimer_types, ONLY: dimer_env_type
16 : USE dimer_utils, ONLY: dimer_thrs, &
17 : rotate_dimer
18 : USE global_types, ONLY: global_environment_type
19 : USE gopt_f_types, ONLY: gopt_f_type
20 : USE gopt_param_types, ONLY: gopt_param_type
21 : USE input_constants, ONLY: default_cell_method_id, &
22 : default_minimization_method_id, &
23 : default_shellcore_method_id, &
24 : default_ts_method_id, &
25 : ls_2pnt, &
26 : ls_fit, &
27 : ls_gold
28 : USE kinds, ONLY: dp
29 : USE mathconstants, ONLY: pi
30 : USE memory_utilities, ONLY: reallocate
31 : USE message_passing, ONLY: mp_para_env_type
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 : PRIVATE
36 :
37 : #:include "gopt_f77_methods.fypp"
38 :
39 : PUBLIC :: cg_linmin, get_conjugate_direction
40 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_utils'
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief Main driver for line minimization routines for CG
47 : !> \param gopt_env ...
48 : !> \param xvec ...
49 : !> \param xi ...
50 : !> \param g ...
51 : !> \param opt_energy ...
52 : !> \param output_unit ...
53 : !> \param gopt_param ...
54 : !> \param globenv ...
55 : !> \par History
56 : !> 10.2005 created [tlaino]
57 : !> \author Teodoro Laino
58 : ! **************************************************************************************************
59 1896 : RECURSIVE SUBROUTINE cg_linmin(gopt_env, xvec, xi, g, opt_energy, output_unit, gopt_param, &
60 : globenv)
61 :
62 : TYPE(gopt_f_type), POINTER :: gopt_env
63 : REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi, g
64 : REAL(KIND=dp), INTENT(INOUT) :: opt_energy
65 : INTEGER :: output_unit
66 : TYPE(gopt_param_type), POINTER :: gopt_param
67 : TYPE(global_environment_type), POINTER :: globenv
68 :
69 : CHARACTER(len=*), PARAMETER :: routineN = 'cg_linmin'
70 :
71 : INTEGER :: handle
72 :
73 1896 : CALL timeset(routineN, handle)
74 1896 : gopt_env%do_line_search = .TRUE.
75 2200 : SELECT CASE (gopt_env%type_id)
76 : CASE (default_minimization_method_id)
77 1328 : SELECT CASE (gopt_param%cg_ls%type_id)
78 : CASE (ls_2pnt)
79 : CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=gopt_param%cg_ls%grad_only, &
80 170 : output_unit=output_unit)
81 : CASE (ls_fit)
82 : CALL linmin_fit(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brack_limit, &
83 54 : gopt_param%cg_ls%initial_step, output_unit, gopt_param, globenv)
84 : CASE (ls_gold)
85 : CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
86 : gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
87 80 : gopt_param%cg_ls%initial_step, output_unit, globenv)
88 : CASE DEFAULT
89 304 : CPABORT("Line Search type not yet implemented in CG.")
90 : END SELECT
91 : CASE (default_ts_method_id)
92 2276 : SELECT CASE (gopt_param%cg_ls%type_id)
93 : CASE (ls_2pnt)
94 854 : IF (gopt_env%dimer_rotation) THEN
95 726 : CALL rotmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy)
96 : ELSE
97 : CALL tslmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy, gopt_param, &
98 128 : output_unit)
99 : END IF
100 : CASE DEFAULT
101 854 : CPABORT("Line Search type not yet implemented in CG for TS search.")
102 : END SELECT
103 : CASE (default_cell_method_id)
104 956 : SELECT CASE (gopt_param%cg_ls%type_id)
105 : CASE (ls_2pnt)
106 : CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.TRUE., &
107 218 : output_unit=output_unit)
108 : CASE (ls_fit)
109 : CALL linmin_fit(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brack_limit, &
110 0 : gopt_param%cg_ls%initial_step, output_unit, gopt_param, globenv)
111 : CASE (ls_gold)
112 : CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
113 : gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
114 350 : gopt_param%cg_ls%initial_step, output_unit, globenv)
115 : CASE DEFAULT
116 568 : CPABORT("Line Search type not yet implemented in CG for cell optimization.")
117 : END SELECT
118 : CASE (default_shellcore_method_id)
119 1896 : SELECT CASE (gopt_param%cg_ls%type_id)
120 : CASE (ls_2pnt)
121 : CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.TRUE., &
122 170 : output_unit=output_unit)
123 : CASE DEFAULT
124 170 : CPABORT("Line Search type not yet implemented in CG for shellcore optimization.")
125 : END SELECT
126 :
127 : END SELECT
128 1896 : gopt_env%do_line_search = .FALSE.
129 1896 : CALL timestop(handle)
130 :
131 1896 : END SUBROUTINE cg_linmin
132 :
133 : ! **************************************************************************************************
134 : !> \brief Line search subroutine based on 2 points (using gradients and energies
135 : !> or only gradients)
136 : !> \param gopt_env ...
137 : !> \param x0 ...
138 : !> \param ls_vec ...
139 : !> \param g ...
140 : !> \param opt_energy ...
141 : !> \param gopt_param ...
142 : !> \param use_only_grad ...
143 : !> \param output_unit ...
144 : !> \author Teodoro Laino - created [tlaino] - 03.2008
145 : ! **************************************************************************************************
146 558 : RECURSIVE SUBROUTINE linmin_2pnt(gopt_env, x0, ls_vec, g, opt_energy, gopt_param, use_only_grad, &
147 : output_unit)
148 : TYPE(gopt_f_type), POINTER :: gopt_env
149 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0, ls_vec, g
150 : REAL(KIND=dp), INTENT(INOUT) :: opt_energy
151 : TYPE(gopt_param_type), POINTER :: gopt_param
152 : LOGICAL, INTENT(IN), OPTIONAL :: use_only_grad
153 : INTEGER, INTENT(IN) :: output_unit
154 :
155 : CHARACTER(len=*), PARAMETER :: routineN = 'linmin_2pnt'
156 :
157 : INTEGER :: handle
158 : LOGICAL :: my_use_only_grad, &
159 : save_consistent_energy_force
160 : REAL(KIND=dp) :: a, b, c, dx, dx_min, dx_min_save, &
161 : dx_thrs, norm_grad1, norm_grad2, &
162 : norm_ls_vec, opt_energy2, x_grad_zero
163 558 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient2, ls_norm
164 :
165 558 : CALL timeset(routineN, handle)
166 298926 : norm_ls_vec = SQRT(DOT_PRODUCT(ls_vec, ls_vec))
167 558 : my_use_only_grad = .FALSE.
168 558 : IF (PRESENT(use_only_grad)) my_use_only_grad = use_only_grad
169 558 : IF (norm_ls_vec /= 0.0_dp) THEN
170 1674 : ALLOCATE (ls_norm(SIZE(ls_vec)))
171 1116 : ALLOCATE (gradient2(SIZE(ls_vec)))
172 597294 : ls_norm = ls_vec/norm_ls_vec
173 558 : dx = norm_ls_vec
174 558 : dx_thrs = gopt_param%cg_ls%max_step
175 :
176 597294 : x0 = x0 + dx*ls_norm
177 : ![NB] don't need consistent energies and forces if using only gradient
178 558 : save_consistent_energy_force = gopt_env%require_consistent_energy_force
179 558 : gopt_env%require_consistent_energy_force = .NOT. my_use_only_grad
180 : CALL cp_eval_at(gopt_env, x0, opt_energy2, gradient2, master=gopt_env%force_env%para_env%mepos, &
181 558 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
182 558 : gopt_env%require_consistent_energy_force = save_consistent_energy_force
183 :
184 298926 : norm_grad1 = -DOT_PRODUCT(g, ls_norm)
185 298926 : norm_grad2 = DOT_PRODUCT(gradient2, ls_norm)
186 558 : IF (my_use_only_grad) THEN
187 : ! a*x+b=y
188 : ! per x=0; b=norm_grad1
189 388 : b = norm_grad1
190 : ! per x=dx; a*dx+b=norm_grad2
191 388 : a = (norm_grad2 - b)/dx
192 388 : x_grad_zero = -b/a
193 388 : dx_min = x_grad_zero
194 : ELSE
195 : ! ax**2+b*x+c=y
196 : ! per x=0 ; c=opt_energy
197 170 : c = opt_energy
198 : ! per x=dx; a*dx**2 + b*dx + c = opt_energy2
199 : ! per x=dx; 2*a*dx + b = norm_grad2
200 : !
201 : ! - a*dx**2 + c = (opt_energy2-norm_grad2*dx)
202 : ! a*dx**2 = c - (opt_energy2-norm_grad2*dx)
203 170 : a = (c - (opt_energy2 - norm_grad2*dx))/dx**2
204 170 : b = norm_grad2 - 2.0_dp*a*dx
205 170 : dx_min = 0.0_dp
206 170 : IF (a /= 0.0_dp) dx_min = -b/(2.0_dp*a)
207 170 : opt_energy = opt_energy2
208 : END IF
209 558 : dx_min_save = dx_min
210 : ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
211 : ! step length
212 558 : IF (ABS(dx_min) > dx_thrs) dx_min = SIGN(1.0_dp, dx_min)*dx_thrs
213 597294 : x0 = x0 + (dx_min - dx)*ls_norm
214 :
215 : ! Print out LS info
216 558 : IF (output_unit > 0) THEN
217 279 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
218 : WRITE (UNIT=output_unit, FMT="(T2,A,T31,A,T78,A)") &
219 279 : "***", "2PNT LINE SEARCH INFO", "***"
220 279 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A)") "***", "***"
221 : WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
222 279 : "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
223 : WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
224 279 : "***", "DX (FITTED )=", dx_min_save, "DX (ACCEPTED )=", dx_min, "***"
225 279 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
226 : END IF
227 558 : DEALLOCATE (ls_norm)
228 558 : DEALLOCATE (gradient2)
229 : ELSE
230 : ! Do Nothing, since.. if the effective force is 0 means that we are already
231 : ! in the saddle point..
232 : END IF
233 558 : CALL timestop(handle)
234 558 : END SUBROUTINE linmin_2pnt
235 :
236 : ! **************************************************************************************************
237 : !> \brief Translational minimization for the Dimer Method - 2pnt LS
238 : !> \param gopt_env ...
239 : !> \param dimer_env ...
240 : !> \param x0 ...
241 : !> \param tls_vec ...
242 : !> \param opt_energy ...
243 : !> \param gopt_param ...
244 : !> \param output_unit ...
245 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
246 : ! **************************************************************************************************
247 128 : SUBROUTINE tslmin_2pnt(gopt_env, dimer_env, x0, tls_vec, opt_energy, gopt_param, output_unit)
248 : TYPE(gopt_f_type), POINTER :: gopt_env
249 : TYPE(dimer_env_type), POINTER :: dimer_env
250 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0, tls_vec
251 : REAL(KIND=dp), INTENT(INOUT) :: opt_energy
252 : TYPE(gopt_param_type), POINTER :: gopt_param
253 : INTEGER, INTENT(IN) :: output_unit
254 :
255 : CHARACTER(len=*), PARAMETER :: routineN = 'tslmin_2pnt'
256 :
257 : INTEGER :: handle
258 : REAL(KIND=dp) :: dx, dx_min, dx_min_acc, dx_min_save, &
259 : dx_thrs, norm_tls_vec, opt_energy2
260 128 : REAL(KIND=dp), DIMENSION(:), POINTER :: tls_norm
261 :
262 128 : CALL timeset(routineN, handle)
263 1886 : norm_tls_vec = SQRT(DOT_PRODUCT(tls_vec, tls_vec))
264 128 : IF (norm_tls_vec /= 0.0_dp) THEN
265 384 : ALLOCATE (tls_norm(SIZE(tls_vec)))
266 :
267 3644 : tls_norm = tls_vec/norm_tls_vec
268 128 : dimer_env%tsl%tls_vec => tls_norm
269 :
270 128 : dx = norm_tls_vec
271 128 : dx_thrs = gopt_param%cg_ls%max_step
272 : ! If curvature is positive let's make the largest step allowed
273 128 : IF (dimer_env%rot%curvature > 0) dx = dx_thrs
274 3644 : x0 = x0 + dx*tls_norm
275 : CALL cp_eval_at(gopt_env, x0, opt_energy2, master=gopt_env%force_env%para_env%mepos, &
276 128 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
277 128 : IF (dimer_env%rot%curvature > 0) THEN
278 30 : dx_min = 0.0_dp
279 30 : dx_min_save = dx
280 30 : dx_min_acc = dx
281 : ELSE
282 : ! First let's try to interpolate the minimum
283 98 : dx_min = -opt_energy/(opt_energy2 - opt_energy)*dx
284 : ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
285 : ! step length
286 98 : dx_min_save = dx_min
287 98 : IF (ABS(dx_min) > dx_thrs) dx_min = SIGN(1.0_dp, dx_min)*dx_thrs
288 98 : dx_min_acc = dx_min
289 98 : dx_min = dx_min - dx
290 : END IF
291 3644 : x0 = x0 + dx_min*tls_norm
292 :
293 : ! Print out LS info
294 128 : IF (output_unit > 0) THEN
295 64 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
296 : WRITE (UNIT=output_unit, FMT="(T2,A,T24,A,T78,A)") &
297 64 : "***", "2PNT TRANSLATIONAL LINE SEARCH INFO", "***"
298 64 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A)") "***", "***"
299 : WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T78,A)") &
300 64 : "***", "LOCAL CURVATURE =", dimer_env%rot%curvature, "***"
301 : WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
302 64 : "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
303 : WRITE (UNIT=output_unit, FMT="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
304 64 : "***", "DX (FITTED )=", dx_min_save, "DX (ACCEPTED )=", dx_min_acc, "***"
305 64 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
306 : END IF
307 :
308 : ! Here we compute the value of the energy in point zero..
309 : CALL cp_eval_at(gopt_env, x0, opt_energy, master=gopt_env%force_env%para_env%mepos, &
310 128 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
311 :
312 128 : DEALLOCATE (tls_norm)
313 : ELSE
314 : ! Do Nothing, since.. if the effective force is 0 means that we are already
315 : ! in the saddle point..
316 : END IF
317 128 : CALL timestop(handle)
318 :
319 128 : END SUBROUTINE tslmin_2pnt
320 :
321 : ! **************************************************************************************************
322 : !> \brief Rotational minimization for the Dimer Method - 2 pnt LS
323 : !> \param gopt_env ...
324 : !> \param dimer_env ...
325 : !> \param x0 ...
326 : !> \param theta ...
327 : !> \param opt_energy ...
328 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
329 : ! **************************************************************************************************
330 726 : SUBROUTINE rotmin_2pnt(gopt_env, dimer_env, x0, theta, opt_energy)
331 : TYPE(gopt_f_type), POINTER :: gopt_env
332 : TYPE(dimer_env_type), POINTER :: dimer_env
333 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0, theta
334 : REAL(KIND=dp), INTENT(INOUT) :: opt_energy
335 :
336 : CHARACTER(len=*), PARAMETER :: routineN = 'rotmin_2pnt'
337 :
338 : INTEGER :: handle
339 : REAL(KIND=dp) :: a0, a1, angle, b1, curvature0, &
340 : curvature1, curvature2, dCdp, f
341 726 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
342 :
343 726 : CALL timeset(routineN, handle)
344 726 : curvature0 = dimer_env%rot%curvature
345 726 : dCdp = dimer_env%rot%dCdp
346 726 : b1 = 0.5_dp*dCdp
347 726 : angle = -0.5_dp*ATAN(dCdp/(2.0_dp*ABS(curvature0)))
348 726 : dimer_env%rot%angle1 = angle
349 12840 : dimer_env%cg_rot%nvec_old = dimer_env%nvec
350 726 : IF (angle > dimer_env%rot%angle_tol) THEN
351 : ! Rotating the dimer of dtheta degrees
352 702 : CALL rotate_dimer(dimer_env%nvec, theta, angle)
353 : ! Re-compute energy, gradients and rotation vector for new R1
354 : CALL cp_eval_at(gopt_env, x0, f, master=gopt_env%force_env%para_env%mepos, &
355 702 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
356 :
357 702 : curvature1 = dimer_env%rot%curvature
358 702 : a1 = (curvature0 - curvature1 + b1*SIN(2.0_dp*angle))/(1.0_dp - COS(2.0_dp*angle))
359 702 : a0 = 2.0_dp*(curvature0 - a1)
360 702 : angle = 0.5_dp*ATAN(b1/a1)
361 702 : curvature2 = a0/2.0_dp + a1*COS(2.0_dp*angle) + b1*SIN(2.0_dp*angle)
362 702 : IF (curvature2 > curvature0) THEN
363 4 : angle = angle + pi/2.0_dp
364 4 : curvature2 = a0/2.0_dp + a1*COS(2.0_dp*angle) + b1*SIN(2.0_dp*angle)
365 : END IF
366 702 : dimer_env%rot%angle2 = angle
367 702 : dimer_env%rot%curvature = curvature2
368 : ! Rotating the dimer the optimized (in plane) vector position
369 12708 : dimer_env%nvec = dimer_env%cg_rot%nvec_old
370 702 : CALL rotate_dimer(dimer_env%nvec, theta, angle)
371 :
372 : ! Evaluate (by interpolation) the norm of the rotational force in the
373 : ! minimum of the rotational search (this is for print-out only)
374 2106 : ALLOCATE (work(SIZE(dimer_env%nvec)))
375 24714 : work = dimer_env%rot%g1
376 : work = SIN(dimer_env%rot%angle1 - dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1 + &
377 : SIN(dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1p + &
378 : (1.0_dp - COS(dimer_env%rot%angle2) - SIN(dimer_env%rot%angle2)*TAN(dimer_env%rot%angle1/2.0_dp))* &
379 24714 : dimer_env%rot%g0
380 24714 : work = -2.0_dp*(work - dimer_env%rot%g0)
381 36720 : work = work - DOT_PRODUCT(work, dimer_env%nvec)*dimer_env%nvec
382 12708 : opt_energy = SQRT(DOT_PRODUCT(work, work))
383 702 : DEALLOCATE (work)
384 : END IF
385 726 : dimer_env%rot%angle2 = angle
386 726 : CALL timestop(handle)
387 :
388 726 : END SUBROUTINE rotmin_2pnt
389 :
390 : ! **************************************************************************************************
391 : !> \brief Line Minimization - Fit
392 : !> \param gopt_env ...
393 : !> \param xvec ...
394 : !> \param xi ...
395 : !> \param opt_energy ...
396 : !> \param brack_limit ...
397 : !> \param step ...
398 : !> \param output_unit ...
399 : !> \param gopt_param ...
400 : !> \param globenv ...
401 : !> \par History
402 : !> 10.2005 created [tlaino]
403 : !> \author Teodoro Laino
404 : !> \note
405 : !> Given as input the vector XVEC and XI, finds the scalar
406 : !> xmin that minimizes the energy XVEC+xmin*XI. Replace step
407 : !> with the optimal value. Enhanced Version
408 : ! **************************************************************************************************
409 54 : SUBROUTINE linmin_fit(gopt_env, xvec, xi, opt_energy, &
410 : brack_limit, step, output_unit, gopt_param, globenv)
411 : TYPE(gopt_f_type), POINTER :: gopt_env
412 : REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi
413 : REAL(KIND=dp) :: opt_energy, brack_limit, step
414 : INTEGER :: output_unit
415 : TYPE(gopt_param_type), POINTER :: gopt_param
416 : TYPE(global_environment_type), POINTER :: globenv
417 :
418 : CHARACTER(len=*), PARAMETER :: routineN = 'linmin_fit'
419 :
420 : INTEGER :: handle, loc_iter, odim
421 : LOGICAL :: should_stop
422 : REAL(KIND=dp) :: ax, bx, fprev, rms_dr, rms_force, scale, &
423 : xmin, xx
424 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
425 54 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hist
426 :
427 54 : CALL timeset(routineN, handle)
428 :
429 54 : NULLIFY (pcom, xicom, hist)
430 54 : rms_dr = gopt_param%rms_dr
431 54 : rms_force = gopt_param%rms_force
432 162 : ALLOCATE (pcom(SIZE(xvec)))
433 108 : ALLOCATE (xicom(SIZE(xvec)))
434 :
435 12342 : pcom = xvec
436 12342 : xicom = xi
437 12342 : xicom = xicom/SQRT(DOT_PRODUCT(xicom, xicom))
438 54 : step = step*0.8_dp ! target a little before the minimum for the first point
439 54 : ax = 0.0_dp
440 54 : xx = step
441 : CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
442 54 : histpoint=hist, globenv=globenv)
443 : !
444 54 : fprev = 0.0_dp
445 234 : opt_energy = MINVAL(hist(:, 2))
446 54 : odim = SIZE(hist, 1)
447 54 : scale = 0.25_dp
448 54 : loc_iter = 0
449 330 : DO WHILE (ABS(hist(odim, 3)) > rms_force*scale .OR. ABS(hist(odim, 1) - hist(odim - 1, 1)) > scale*rms_dr)
450 276 : CALL external_control(should_stop, "LINFIT", globenv=globenv)
451 276 : IF (should_stop) EXIT
452 : !
453 276 : loc_iter = loc_iter + 1
454 276 : fprev = opt_energy
455 276 : xmin = FindMin(hist(:, 1), hist(:, 2), hist(:, 3))
456 276 : CALL reallocate(hist, 1, odim + 1, 1, 3)
457 276 : hist(odim + 1, 1) = xmin
458 276 : hist(odim + 1, 3) = cg_deval1d(gopt_env, xmin, pcom, xicom, opt_energy)
459 276 : hist(odim + 1, 2) = opt_energy
460 276 : odim = SIZE(hist, 1)
461 : END DO
462 : !
463 6198 : xicom = xmin*xicom
464 54 : step = xmin
465 12342 : xvec = xvec + xicom
466 54 : DEALLOCATE (pcom)
467 54 : DEALLOCATE (xicom)
468 54 : DEALLOCATE (hist)
469 54 : IF (output_unit > 0) THEN
470 27 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
471 : WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
472 27 : "***", "FIT LS - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
473 27 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
474 : END IF
475 54 : CALL timestop(handle)
476 :
477 54 : END SUBROUTINE linmin_fit
478 :
479 : ! **************************************************************************************************
480 : !> \brief Line Minimization routine - GOLD
481 : !> \param gopt_env ...
482 : !> \param xvec ...
483 : !> \param xi ...
484 : !> \param opt_energy ...
485 : !> \param brent_tol ...
486 : !> \param brent_max_iter ...
487 : !> \param brack_limit ...
488 : !> \param step ...
489 : !> \param output_unit ...
490 : !> \param globenv ...
491 : !> \par History
492 : !> 10.2005 created [tlaino]
493 : !> \author Teodoro Laino
494 : !> \note
495 : !> Given as input the vector XVEC and XI, finds the scalar
496 : !> xmin that minimizes the energy XVEC+xmin*XI. Replaces XMIN
497 : !> with the optimal value
498 : ! **************************************************************************************************
499 430 : SUBROUTINE linmin_gold(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
500 : brack_limit, step, output_unit, globenv)
501 : TYPE(gopt_f_type), POINTER :: gopt_env
502 : REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi
503 : REAL(KIND=dp) :: opt_energy, brent_tol
504 : INTEGER :: brent_max_iter
505 : REAL(KIND=dp) :: brack_limit, step
506 : INTEGER :: output_unit
507 : TYPE(global_environment_type), POINTER :: globenv
508 :
509 : CHARACTER(len=*), PARAMETER :: routineN = 'linmin_gold'
510 :
511 : INTEGER :: handle
512 : REAL(KIND=dp) :: ax, bx, xmin, xx
513 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
514 :
515 430 : CALL timeset(routineN, handle)
516 :
517 : NULLIFY (pcom, xicom)
518 1290 : ALLOCATE (pcom(SIZE(xvec)))
519 860 : ALLOCATE (xicom(SIZE(xvec)))
520 :
521 122686 : pcom = xvec
522 122686 : xicom = xi
523 122686 : xicom = xicom/SQRT(DOT_PRODUCT(xicom, xicom))
524 430 : step = step*0.8_dp ! target a little before the minimum for the first point
525 430 : ax = 0.0_dp
526 430 : xx = step
527 : CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
528 430 : globenv=globenv)
529 :
530 : opt_energy = cg_dbrent(gopt_env, ax, xx, bx, brent_tol, brent_max_iter, &
531 430 : xmin, pcom, xicom, output_unit, globenv)
532 61558 : xicom = xmin*xicom
533 430 : step = xmin
534 122686 : xvec = xvec + xicom
535 430 : DEALLOCATE (pcom)
536 430 : DEALLOCATE (xicom)
537 430 : CALL timestop(handle)
538 430 : END SUBROUTINE linmin_gold
539 :
540 : ! **************************************************************************************************
541 : !> \brief Routine for initially bracketing a minimum based on the golden search
542 : !> minimum
543 : !> \param gopt_env ...
544 : !> \param ax ...
545 : !> \param bx ...
546 : !> \param cx ...
547 : !> \param pcom ...
548 : !> \param xicom ...
549 : !> \param brack_limit ...
550 : !> \param output_unit ...
551 : !> \param histpoint ...
552 : !> \param globenv ...
553 : !> \par History
554 : !> 10.2005 created [tlaino]
555 : !> \author Teodoro Laino
556 : !> \note
557 : !> Given two distinct initial points ax and bx this routine searches
558 : !> in the downhill direction and returns new points ax, bx, cx that
559 : !> bracket the minimum of the function
560 : ! **************************************************************************************************
561 484 : SUBROUTINE cg_mnbrak(gopt_env, ax, bx, cx, pcom, xicom, brack_limit, output_unit, &
562 : histpoint, globenv)
563 : TYPE(gopt_f_type), POINTER :: gopt_env
564 : REAL(KIND=dp) :: ax, bx, cx
565 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
566 : REAL(KIND=dp) :: brack_limit
567 : INTEGER :: output_unit
568 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: histpoint
569 : TYPE(global_environment_type), POINTER :: globenv
570 :
571 : CHARACTER(len=*), PARAMETER :: routineN = 'cg_mnbrak'
572 :
573 : INTEGER :: handle, loc_iter, odim
574 : LOGICAL :: hist, should_stop
575 : REAL(KIND=dp) :: dum, fa, fb, fc, fu, gold, q, r, u, ulim
576 :
577 484 : CALL timeset(routineN, handle)
578 484 : hist = PRESENT(histpoint)
579 484 : IF (hist) THEN
580 54 : CPASSERT(.NOT. ASSOCIATED(histpoint))
581 54 : ALLOCATE (histpoint(3, 3))
582 : END IF
583 484 : gold = (1.0_dp + SQRT(5.0_dp))/2.0_dp
584 : IF (hist) THEN
585 54 : histpoint(1, 1) = ax
586 54 : histpoint(1, 3) = cg_deval1d(gopt_env, ax, pcom, xicom, fa)
587 54 : histpoint(1, 2) = fa
588 54 : histpoint(2, 1) = bx
589 54 : histpoint(2, 3) = cg_deval1d(gopt_env, bx, pcom, xicom, fb)
590 54 : histpoint(2, 2) = fb
591 : ELSE
592 430 : fa = cg_eval1d(gopt_env, ax, pcom, xicom)
593 430 : fb = cg_eval1d(gopt_env, bx, pcom, xicom)
594 : END IF
595 484 : IF (fb > fa) THEN
596 138 : dum = ax
597 138 : ax = bx
598 138 : bx = dum
599 138 : dum = fb
600 138 : fb = fa
601 138 : fa = dum
602 : END IF
603 484 : cx = bx + gold*(bx - ax)
604 484 : IF (hist) THEN
605 54 : histpoint(3, 1) = cx
606 54 : histpoint(3, 3) = cg_deval1d(gopt_env, cx, pcom, xicom, fc)
607 54 : histpoint(3, 2) = fc
608 : ELSE
609 430 : fc = cg_eval1d(gopt_env, cx, pcom, xicom)
610 : END IF
611 484 : loc_iter = 3
612 588 : DO WHILE (fb >= fc)
613 146 : CALL external_control(should_stop, "MNBRACK", globenv=globenv)
614 146 : IF (should_stop) EXIT
615 : !
616 146 : r = (bx - ax)*(fb - fc)
617 146 : q = (bx - cx)*(fb - fa)
618 146 : u = bx - ((bx - cx)*q - (bx - ax)*r)/(2.0_dp*SIGN(MAX(ABS(q - r), TINY(0.0_dp)), q - r))
619 146 : ulim = bx + brack_limit*(cx - bx)
620 146 : IF ((bx - u)*(u - cx) > 0.0_dp) THEN
621 46 : IF (hist) THEN
622 10 : odim = SIZE(histpoint, 1)
623 10 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
624 10 : histpoint(odim + 1, 1) = u
625 10 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
626 10 : histpoint(odim + 1, 2) = fu
627 : ELSE
628 36 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
629 : END IF
630 46 : loc_iter = loc_iter + 1
631 46 : IF (fu < fc) THEN
632 42 : ax = bx
633 : fa = fb
634 42 : bx = u
635 : fb = fu
636 42 : EXIT
637 4 : ELSE IF (fu > fb) THEN
638 0 : cx = u
639 : fc = fu
640 0 : EXIT
641 : END IF
642 4 : u = cx + gold*(cx - bx)
643 4 : IF (hist) THEN
644 0 : odim = SIZE(histpoint, 1)
645 0 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
646 0 : histpoint(odim + 1, 1) = u
647 0 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
648 0 : histpoint(odim + 1, 2) = fu
649 : ELSE
650 4 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
651 : END IF
652 4 : loc_iter = loc_iter + 1
653 100 : ELSE IF ((cx - u)*(u - ulim) > 0.) THEN
654 100 : IF (hist) THEN
655 4 : odim = SIZE(histpoint, 1)
656 4 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
657 4 : histpoint(odim + 1, 1) = u
658 4 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
659 4 : histpoint(odim + 1, 2) = fu
660 : ELSE
661 96 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
662 : END IF
663 100 : loc_iter = loc_iter + 1
664 100 : IF (fu < fc) THEN
665 98 : bx = cx
666 98 : cx = u
667 98 : u = cx + gold*(cx - bx)
668 98 : fb = fc
669 98 : fc = fu
670 98 : IF (hist) THEN
671 4 : odim = SIZE(histpoint, 1)
672 4 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
673 4 : histpoint(odim + 1, 1) = u
674 4 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
675 4 : histpoint(odim + 1, 2) = fu
676 : ELSE
677 94 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
678 : END IF
679 98 : loc_iter = loc_iter + 1
680 : END IF
681 0 : ELSE IF ((u - ulim)*(ulim - cx) >= 0.) THEN
682 0 : u = ulim
683 0 : IF (hist) THEN
684 0 : odim = SIZE(histpoint, 1)
685 0 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
686 0 : histpoint(odim + 1, 1) = u
687 0 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
688 0 : histpoint(odim + 1, 2) = fu
689 : ELSE
690 0 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
691 : END IF
692 0 : loc_iter = loc_iter + 1
693 : ELSE
694 0 : u = cx + gold*(cx - bx)
695 0 : IF (hist) THEN
696 0 : odim = SIZE(histpoint, 1)
697 0 : CALL reallocate(histpoint, 1, odim + 1, 1, 3)
698 0 : histpoint(odim + 1, 1) = u
699 0 : histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
700 0 : histpoint(odim + 1, 2) = fu
701 : ELSE
702 0 : fu = cg_eval1d(gopt_env, u, pcom, xicom)
703 : END IF
704 0 : loc_iter = loc_iter + 1
705 : END IF
706 104 : ax = bx
707 104 : bx = cx
708 104 : cx = u
709 104 : fa = fb
710 104 : fb = fc
711 104 : fc = fu
712 : END DO
713 484 : IF (output_unit > 0) THEN
714 242 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
715 : WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
716 242 : "***", "MNBRACK - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
717 242 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
718 : END IF
719 484 : CALL timestop(handle)
720 484 : END SUBROUTINE cg_mnbrak
721 :
722 : ! **************************************************************************************************
723 : !> \brief Routine implementing the Brent Method
724 : !> Brent,R.P. Algorithm for Minimization without Derivatives, Chapt.5
725 : !> 1973
726 : !> Extension in the use of derivatives
727 : !> \param gopt_env ...
728 : !> \param ax ...
729 : !> \param bx ...
730 : !> \param cx ...
731 : !> \param tol ...
732 : !> \param itmax ...
733 : !> \param xmin ...
734 : !> \param pcom ...
735 : !> \param xicom ...
736 : !> \param output_unit ...
737 : !> \param globenv ...
738 : !> \return ...
739 : !> \par History
740 : !> 10.2005 created [tlaino]
741 : !> \author Teodoro Laino
742 : !> \note
743 : !> Given a bracketing triplet of abscissas ax, bx, cx (such that bx
744 : !> is between ax and cx and energy of bx is less than energy of ax and cx),
745 : !> this routine isolates the minimum to a precision of about tol using
746 : !> Brent method. This routine implements the extension of the Brent Method
747 : !> using derivatives
748 : ! **************************************************************************************************
749 430 : FUNCTION cg_dbrent(gopt_env, ax, bx, cx, tol, itmax, xmin, pcom, xicom, output_unit, &
750 : globenv) RESULT(dbrent)
751 : TYPE(gopt_f_type), POINTER :: gopt_env
752 : REAL(KIND=dp) :: ax, bx, cx, tol
753 : INTEGER :: itmax
754 : REAL(KIND=dp) :: xmin
755 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
756 : INTEGER :: output_unit
757 : TYPE(global_environment_type), POINTER :: globenv
758 : REAL(KIND=dp) :: dbrent
759 :
760 : CHARACTER(len=*), PARAMETER :: routineN = 'cg_dbrent'
761 : REAL(KIND=dp), PARAMETER :: zeps = 1.0E-8_dp
762 :
763 : INTEGER :: handle, iter, loc_iter
764 : LOGICAL :: ok1, ok2, should_stop, skip0, skip1
765 : REAL(KIND=dp) :: a, b, d, d1, d2, du, dv, dw, dx, e, fu, &
766 : fv, fw, fx, olde, tol1, tol2, u, u1, &
767 : u2, v, w, x, xm
768 :
769 430 : CALL timeset(routineN, handle)
770 430 : a = MIN(ax, cx)
771 430 : b = MAX(ax, cx)
772 430 : v = bx; w = v; x = v
773 430 : e = 0.0_dp
774 430 : dx = cg_deval1d(gopt_env, x, pcom, xicom, fx)
775 430 : fv = fx
776 430 : fw = fx
777 430 : dv = dx
778 430 : dw = dx
779 430 : loc_iter = 1
780 2086 : DO iter = 1, itmax
781 2086 : CALL external_control(should_stop, "BRENT", globenv=globenv)
782 2086 : IF (should_stop) EXIT
783 : !
784 2086 : xm = 0.5_dp*(a + b)
785 2086 : tol1 = tol*ABS(x) + zeps
786 2086 : tol2 = 2.0_dp*tol1
787 2086 : skip0 = .FALSE.
788 2086 : skip1 = .FALSE.
789 2086 : IF (ABS(x - xm) <= (tol2 - 0.5_dp*(b - a))) EXIT
790 2018 : IF (ABS(e) > tol1) THEN
791 1538 : d1 = 2.0_dp*(b - a)
792 1538 : d2 = d1
793 1538 : IF (dw /= dx) d1 = (w - x)*dx/(dx - dw)
794 1538 : IF (dv /= dx) d2 = (v - x)*dx/(dx - dv)
795 1538 : u1 = x + d1
796 1538 : u2 = x + d2
797 1538 : ok1 = ((a - u1)*(u1 - b) > 0.0_dp) .AND. (dx*d1 <= 0.0_dp)
798 1538 : ok2 = ((a - u2)*(u2 - b) > 0.0_dp) .AND. (dx*d2 <= 0.0_dp)
799 1732 : olde = e
800 1732 : e = d
801 1034 : IF (.NOT. (ok1 .OR. ok2)) THEN
802 : skip0 = .TRUE.
803 698 : ELSE IF (ok1 .AND. ok2) THEN
804 500 : IF (ABS(d1) < ABS(d2)) THEN
805 : d = d1
806 : ELSE
807 : d = d2
808 : END IF
809 198 : ELSE IF (ok1) THEN
810 : d = d1
811 : ELSE
812 : d = d2
813 : END IF
814 : IF (.NOT. skip0) THEN
815 698 : IF (ABS(d) > ABS(0.5_dp*olde)) skip0 = .TRUE.
816 : IF (.NOT. skip0) THEN
817 672 : u = x + d
818 672 : IF ((u - a) < tol2 .OR. (b - u) < tol2) d = SIGN(tol1, xm - x)
819 : skip1 = .TRUE.
820 : END IF
821 : END IF
822 : END IF
823 : IF (.NOT. skip1) THEN
824 1346 : IF (dx >= 0.0_dp) THEN
825 148 : e = a - x
826 : ELSE
827 1198 : e = b - x
828 : END IF
829 1346 : d = 0.5_dp*e
830 : END IF
831 2018 : IF (ABS(d) >= tol1) THEN
832 1574 : u = x + d
833 1574 : du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
834 1574 : loc_iter = loc_iter + 1
835 : ELSE
836 444 : u = x + SIGN(tol1, d)
837 444 : du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
838 444 : loc_iter = loc_iter + 1
839 444 : IF (fu > fx) EXIT
840 : END IF
841 4172 : IF (fu <= fx) THEN
842 626 : IF (u >= x) THEN
843 : a = x
844 : ELSE
845 290 : b = x
846 : END IF
847 626 : v = w; fv = fw; dv = dw; w = x
848 626 : fw = fx; dw = dx; x = u; fx = fu; dx = du
849 : ELSE
850 1030 : IF (u < x) THEN
851 : a = u
852 : ELSE
853 918 : b = u
854 : END IF
855 1030 : IF (fu <= fw .OR. w == x) THEN
856 : v = w; fv = fw; dv = dw
857 : w = u; fw = fu; dw = du
858 140 : ELSE IF (fu <= fv .OR. v == x .OR. v == w) THEN
859 140 : v = u
860 140 : fv = fu
861 140 : dv = du
862 : END IF
863 : END IF
864 : END DO
865 430 : IF (output_unit > 0) THEN
866 215 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
867 : WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,I7,T78,A)") &
868 215 : "***", "BRENT - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
869 215 : IF (iter == itmax + 1) &
870 : WRITE (UNIT=output_unit, FMT="(T2,A,T22,A,T78,A)") &
871 0 : "***", "BRENT - NUMBER OF ITERATIONS EXCEEDED ", "***"
872 215 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
873 : END IF
874 430 : CPASSERT(iter /= itmax + 1)
875 430 : xmin = x
876 430 : dbrent = fx
877 430 : CALL timestop(handle)
878 :
879 430 : END FUNCTION cg_dbrent
880 :
881 : ! **************************************************************************************************
882 : !> \brief Evaluates energy in one dimensional space defined by the point
883 : !> pcom and with direction xicom, position x
884 : !> \param gopt_env ...
885 : !> \param x ...
886 : !> \param pcom ...
887 : !> \param xicom ...
888 : !> \return ...
889 : !> \par History
890 : !> 10.2005 created [tlaino]
891 : !> \author Teodoro Laino
892 : ! **************************************************************************************************
893 1520 : FUNCTION cg_eval1d(gopt_env, x, pcom, xicom) RESULT(my_val)
894 : TYPE(gopt_f_type), POINTER :: gopt_env
895 : REAL(KIND=dp) :: x
896 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
897 : REAL(KIND=dp) :: my_val
898 :
899 : CHARACTER(len=*), PARAMETER :: routineN = 'cg_eval1d'
900 :
901 : INTEGER :: handle
902 1520 : REAL(KIND=dp), DIMENSION(:), POINTER :: xvec
903 :
904 1520 : CALL timeset(routineN, handle)
905 :
906 4560 : ALLOCATE (xvec(SIZE(pcom)))
907 441280 : xvec = pcom + x*xicom
908 : CALL cp_eval_at(gopt_env, xvec, my_val, master=gopt_env%force_env%para_env%mepos, &
909 1520 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
910 1520 : DEALLOCATE (xvec)
911 :
912 1520 : CALL timestop(handle)
913 :
914 1520 : END FUNCTION cg_eval1d
915 :
916 : ! **************************************************************************************************
917 : !> \brief Evaluates derivatives in one dimensional space defined by the point
918 : !> pcom and with direction xicom, position x
919 : !> \param gopt_env ...
920 : !> \param x ...
921 : !> \param pcom ...
922 : !> \param xicom ...
923 : !> \param fval ...
924 : !> \return ...
925 : !> \par History
926 : !> 10.2005 created [tlaino]
927 : !> \author Teodoro Laino
928 : ! **************************************************************************************************
929 2904 : FUNCTION cg_deval1d(gopt_env, x, pcom, xicom, fval) RESULT(my_val)
930 : TYPE(gopt_f_type), POINTER :: gopt_env
931 : REAL(KIND=dp) :: x
932 : REAL(KIND=dp), DIMENSION(:), POINTER :: pcom, xicom
933 : REAL(KIND=dp) :: fval, my_val
934 :
935 : CHARACTER(len=*), PARAMETER :: routineN = 'cg_deval1d'
936 :
937 : INTEGER :: handle
938 : REAL(KIND=dp) :: energy
939 2904 : REAL(KIND=dp), DIMENSION(:), POINTER :: grad, xvec
940 :
941 2904 : CALL timeset(routineN, handle)
942 :
943 8712 : ALLOCATE (xvec(SIZE(pcom)))
944 5808 : ALLOCATE (grad(SIZE(pcom)))
945 719280 : xvec = pcom + x*xicom
946 : CALL cp_eval_at(gopt_env, xvec, energy, grad, master=gopt_env%force_env%para_env%mepos, &
947 2904 : final_evaluation=.FALSE., para_env=gopt_env%force_env%para_env)
948 359640 : my_val = DOT_PRODUCT(grad, xicom)
949 2904 : fval = energy
950 2904 : DEALLOCATE (xvec)
951 2904 : DEALLOCATE (grad)
952 2904 : CALL timestop(handle)
953 :
954 2904 : END FUNCTION cg_deval1d
955 :
956 : ! **************************************************************************************************
957 : !> \brief Find the minimum of a parabolic function obtained with a least square fit
958 : !> \param x ...
959 : !> \param y ...
960 : !> \param dy ...
961 : !> \return ...
962 : !> \par History
963 : !> 10.2005 created [fawzi]
964 : !> \author Fawzi Mohamed
965 : ! **************************************************************************************************
966 276 : FUNCTION FindMin(x, y, dy) RESULT(res)
967 : REAL(kind=dp), DIMENSION(:) :: x, y, dy
968 : REAL(kind=dp) :: res
969 :
970 : INTEGER :: i, info, iwork(8*3), lwork, min_pos, np
971 : REAL(kind=dp) :: diag(3), res1(3), res2(3), res3(3), &
972 : spread, sum_x, sum_xx, tmpw(1), &
973 : vt(3, 3)
974 276 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
975 552 : REAL(kind=dp), DIMENSION(2*SIZE(x), 3) :: f
976 552 : REAL(kind=dp), DIMENSION(2*SIZE(x)) :: b, w
977 552 : REAL(kind=dp) :: u(2*SIZE(x), 3)
978 :
979 276 : np = SIZE(x)
980 276 : CPASSERT(np > 1)
981 276 : sum_x = 0._dp
982 276 : sum_xx = 0._dp
983 276 : min_pos = 1
984 2178 : DO i = 1, np
985 1902 : sum_xx = sum_xx + x(i)**2
986 1902 : sum_x = sum_x + x(i)
987 2178 : IF (y(min_pos) > y(i)) min_pos = i
988 : END DO
989 : spread = SQRT(sum_xx/REAL(np, dp) - (sum_x/REAL(np, dp))**2)
990 2178 : DO i = 1, np
991 1902 : w(i) = EXP(-(REAL(np - i, dp))**2/(REAL(2*9, dp)))
992 2178 : w(i + np) = 2._dp*w(i)
993 : END DO
994 2178 : DO i = 1, np
995 1902 : f(i, 1) = w(i)
996 1902 : f(i, 2) = x(i)*w(i)
997 1902 : f(i, 3) = x(i)**2*w(i)
998 1902 : f(i + np, 1) = 0
999 1902 : f(i + np, 2) = w(i + np)
1000 2178 : f(i + np, 3) = 2*x(i)*w(i + np)
1001 : END DO
1002 2178 : DO i = 1, np
1003 1902 : b(i) = y(i)*w(i)
1004 2178 : b(i + np) = dy(i)*w(i + np)
1005 : END DO
1006 276 : lwork = -1
1007 : CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), tmpw, lwork, &
1008 276 : iwork, info)
1009 276 : lwork = CEILING(tmpw(1))
1010 828 : ALLOCATE (work(lwork))
1011 : CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), work, lwork, &
1012 276 : iwork, info)
1013 276 : DEALLOCATE (work)
1014 276 : CALL dgemv('T', SIZE(u, 1), SIZE(u, 2), 1._dp, u, SIZE(u, 1), b, 1, 0._dp, res1, 1)
1015 1104 : DO i = 1, 3
1016 1104 : res2(i) = res1(i)/diag(i)
1017 : END DO
1018 276 : CALL dgemv('T', 3, 3, 1._dp, vt, SIZE(vt, 1), res2, 1, 0._dp, res3, 1)
1019 276 : res = -0.5*res3(2)/res3(3)
1020 276 : END FUNCTION FindMin
1021 :
1022 : ! **************************************************************************************************
1023 : !> \brief Computes the Conjugate direction for the next search
1024 : !> \param gopt_env ...
1025 : !> \param Fletcher_Reeves ...
1026 : !> \param g contains the theta of the previous step.. (norm 1.0 vector)
1027 : !> \param xi contains the -theta of the present step.. (norm 1.0 vector)
1028 : !> \param h contains the search direction of the previous step (must be orthogonal
1029 : !> to nvec of the previous step (nvec_old))
1030 : !> \par Info for DIMER method
1031 : !> \par History
1032 : !> 10.2005 created [tlaino]
1033 : !> \author Teodoro Laino
1034 : ! **************************************************************************************************
1035 1634 : SUBROUTINE get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
1036 : TYPE(gopt_f_type), POINTER :: gopt_env
1037 : LOGICAL, INTENT(IN) :: Fletcher_Reeves
1038 : REAL(KIND=dp), DIMENSION(:), POINTER :: g, xi, h
1039 :
1040 : CHARACTER(len=*), PARAMETER :: routineN = 'get_conjugate_direction'
1041 :
1042 : INTEGER :: handle
1043 : LOGICAL :: check
1044 : REAL(KIND=dp) :: dgg, gam, gg, norm, norm_h
1045 : TYPE(dimer_env_type), POINTER :: dimer_env
1046 :
1047 1634 : CALL timeset(routineN, handle)
1048 1634 : NULLIFY (dimer_env)
1049 1634 : IF (.NOT. gopt_env%dimer_rotation) THEN
1050 298388 : gg = DOT_PRODUCT(g, g)
1051 1040 : IF (Fletcher_Reeves) THEN
1052 0 : dgg = DOT_PRODUCT(xi, xi)
1053 : ELSE
1054 298388 : dgg = DOT_PRODUCT((xi + g), xi)
1055 : END IF
1056 1040 : gam = dgg/gg
1057 595736 : g = h
1058 595736 : h = -xi + gam*h
1059 : ELSE
1060 594 : dimer_env => gopt_env%dimer_env
1061 10932 : check = ABS(DOT_PRODUCT(g, g) - 1.0_dp) < MAX(1.0E-9_dp, dimer_thrs)
1062 594 : CPASSERT(check)
1063 :
1064 10932 : check = ABS(DOT_PRODUCT(xi, xi) - 1.0_dp) < MAX(1.0E-9_dp, dimer_thrs)
1065 594 : CPASSERT(check)
1066 :
1067 10932 : check = ABS(DOT_PRODUCT(h, dimer_env%cg_rot%nvec_old)) < MAX(1.0E-9_dp, dimer_thrs)
1068 594 : CPASSERT(check)
1069 594 : gg = dimer_env%cg_rot%norm_theta_old**2
1070 594 : IF (Fletcher_Reeves) THEN
1071 0 : dgg = dimer_env%cg_rot%norm_theta**2
1072 : ELSE
1073 594 : norm = dimer_env%cg_rot%norm_theta*dimer_env%cg_rot%norm_theta_old
1074 10932 : dgg = dimer_env%cg_rot%norm_theta**2 + DOT_PRODUCT(g, xi)*norm
1075 : END IF
1076 : ! Compute Theta** and store it in nvec_old
1077 594 : CALL rotate_dimer(dimer_env%cg_rot%nvec_old, g, dimer_env%rot%angle2 + pi/2.0_dp)
1078 594 : gam = dgg/gg
1079 21270 : g = h
1080 21270 : h = -xi*dimer_env%cg_rot%norm_theta + gam*dimer_env%cg_rot%norm_h*dimer_env%cg_rot%nvec_old
1081 31608 : h = h - DOT_PRODUCT(h, dimer_env%nvec)*dimer_env%nvec
1082 10932 : norm_h = SQRT(DOT_PRODUCT(h, h))
1083 594 : IF (norm_h < EPSILON(0.0_dp)) THEN
1084 0 : h = 0.0_dp
1085 : ELSE
1086 10932 : h = h/norm_h
1087 : END IF
1088 594 : dimer_env%cg_rot%norm_h = norm_h
1089 : END IF
1090 1634 : CALL timestop(handle)
1091 :
1092 1634 : END SUBROUTINE get_conjugate_direction
1093 :
1094 : END MODULE cg_utils
|
