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