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