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