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