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