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