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