1! 2! Copyright (C) 2003-2010 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!---------------------------------------------------------------------------- 9MODULE bfgs_module 10 !---------------------------------------------------------------------------- 11 ! 12 ! ... Ionic relaxation through the Newton-Raphson optimization scheme 13 ! ... based on the Broyden-Fletcher-Goldfarb-Shanno algorithm for the 14 ! ... estimate of the inverse Hessian matrix. 15 ! ... The ionic relaxation is performed converting cartesian (and cell) 16 ! ... positions into internal coordinates. 17 ! ... The algorithm uses a "trust radius" line search based on Wolfe 18 ! ... conditions. Steps are rejected until the first Wolfe condition 19 ! ... (sufficient energy decrease) is satisfied. Updated step length 20 ! ... is estimated from quadratic interpolation. 21 ! ... When the step is accepted inverse hessian is updated according to 22 ! ... BFGS scheme and a new search direction is obtained from NR or GDIIS 23 ! ... method. The corresponding step length is limited by trust_radius_max 24 ! ... and can't be larger than the previous step multiplied by a certain 25 ! ... factor determined by Wolfe and other convergence conditions. 26 ! 27 ! ... Originally written ( 5/12/2003 ) and maintained ( 2003-2007 ) by 28 ! ... Carlo Sbraccia 29 ! ... Modified for variable-cell-shape relaxation ( 2007-2008 ) by 30 ! ... Javier Antonio Montoya, Lorenzo Paulatto and Stefano de Gironcoli 31 ! ... Re-analyzed by Stefano de Gironcoli ( 2010 ) 32 ! 33 ! ... references : 34 ! 35 ! ... 1) Roger Fletcher, Practical Methods of Optimization, John Wiley and 36 ! ... Sons, Chichester, 2nd edn, 1987. 37 ! ... 2) Salomon R. Billeter, Alexander J. Turner, Walter Thiel, 38 ! ... Phys. Chem. Chem. Phys. 2, 2177 (2000). 39 ! ... 3) Salomon R. Billeter, Alessandro Curioni, Wanda Andreoni, 40 ! ... Comput. Mat. Science 27, 437, (2003). 41 ! ... 4) Ren Weiqing, PhD Thesis: Numerical Methods for the Study of Energy 42 ! ... Landscapes and Rare Events. 43 ! 44 ! 45 USE kinds, ONLY : DP 46 USE io_files, ONLY : iunbfgs, prefix 47 USE constants, ONLY : eps8, eps16 48 USE cell_base, ONLY : iforceh ! FIXME: should be passed as argument 49 ! 50 USE basic_algebra_routines 51 USE matrix_inversion 52 ! 53 IMPLICIT NONE 54 ! 55 PRIVATE 56 ! 57 ! ... public methods 58 ! 59 PUBLIC :: bfgs, terminate_bfgs, bfgs_get_n_iter 60 ! 61 ! ... public variables 62 ! 63 PUBLIC :: bfgs_ndim, & 64 trust_radius_ini, trust_radius_min, trust_radius_max, & 65 w_1, w_2 66 ! 67 ! ... global module variables 68 ! 69 SAVE 70 ! 71 CHARACTER (len=8) :: fname="energy" ! name of the function to be minimized 72 ! 73 REAL(DP), ALLOCATABLE :: & 74 pos(:), &! positions + cell 75 grad(:), &! gradients + cell_force 76 pos_p(:), &! positions at the previous accepted iteration 77 grad_p(:), &! gradients at the previous accepted iteration 78 inv_hess(:,:), &! inverse hessian matrix (updated using BFGS formula) 79 metric(:,:), & 80 h_block(:,:), & 81 hinv_block(:,:), & 82 step(:), &! the (new) search direction (normalized NR step) 83 step_old(:), &! the previous search direction (normalized NR step) 84 pos_old(:,:), &! list of m old positions - used only by gdiis 85 grad_old(:,:), &! list of m old gradients - used only by gdiis 86 pos_best(:) ! best extrapolated positions - used only by gdiis 87 REAL(DP) :: & 88 nr_step_length, &! length of (new) Newton-Raphson step 89 nr_step_length_old,&! length of previous Newton-Raphson step 90 trust_radius, &! new displacement along the search direction 91 trust_radius_old, &! old displacement along the search direction 92 energy_p ! energy at previous accepted iteration 93 INTEGER :: & 94 scf_iter, &! number of scf iterations 95 bfgs_iter, &! number of bfgs iterations 96 gdiis_iter, &! number of gdiis iterations 97 tr_min_hit = 0 ! set to 1 if the trust_radius has already been 98 ! set to the minimum value at the previous step 99 ! set to 2 if trust_radius is reset again: exit 100 LOGICAL :: & 101 conv_bfgs ! .TRUE. when bfgs convergence has been achieved 102 ! 103 ! ... default values for the following variables are set in 104 ! ... Modules/read_namelist.f90 (SUBROUTINE ions_defaults) 105 ! 106 ! ... Note that trust_radius_max, trust_radius_min, trust_radius_ini, 107 ! ... w_1, w_2, bfgs_ndim have a default value, but can also be assigned 108 ! ... in the input. 109 ! 110 INTEGER :: & 111 bfgs_ndim ! dimension of the subspace for GDIIS 112 ! fixed to 1 for standard BFGS algorithm 113 REAL(DP) :: & 114 trust_radius_ini, &! suggested initial displacement 115 trust_radius_min, &! minimum allowed displacement 116 trust_radius_max ! maximum allowed displacement 117 118 REAL(DP) :: &! parameters for Wolfe conditions 119 w_1, &! 1st Wolfe condition: sufficient energy decrease 120 w_2 ! 2nd Wolfe condition: sufficient gradient decrease 121 ! 122CONTAINS 123 ! 124 !------------------------------------------------------------------------ 125 SUBROUTINE bfgs( pos_in, h, energy, grad_in, fcell, fixion, scratch, stdout,& 126 energy_thr, grad_thr, cell_thr, energy_error, grad_error, & 127 cell_error, lmovecell, step_accepted, stop_bfgs, istep ) 128 !------------------------------------------------------------------------ 129 ! 130 ! ... list of input/output arguments : 131 ! 132 ! pos : vector containing 3N coordinates of the system ( x ) 133 ! energy : energy of the system ( V(x) ) 134 ! grad : vector containing 3N components of grad( V(x) ) 135 ! fixion : vector used to freeze a deg. of freedom 136 ! scratch : scratch directory 137 ! stdout : unit for standard output 138 ! energy_thr : treshold on energy difference for BFGS convergence 139 ! grad_thr : treshold on grad difference for BFGS convergence 140 ! the largest component of grad( V(x) ) is considered 141 ! energy_error : energy difference | V(x_i) - V(x_i-1) | 142 ! grad_error : the largest component of 143 ! | grad(V(x_i)) - grad(V(x_i-1)) | 144 ! cell_error : the largest component of: omega*(stress-press*I) 145 ! step_accepted : .TRUE. if a new BFGS step is done 146 ! stop_bfgs : .TRUE. if BFGS convergence has been achieved 147 ! 148 IMPLICIT NONE 149 ! 150 REAL(DP), INTENT(INOUT) :: pos_in(:) 151 REAL(DP), INTENT(INOUT) :: h(3,3) 152 REAL(DP), INTENT(INOUT) :: energy 153 REAL(DP), INTENT(INOUT) :: grad_in(:) 154 REAL(DP), INTENT(INOUT) :: fcell(3,3) 155 INTEGER, INTENT(IN) :: fixion(:) 156 CHARACTER(LEN=*), INTENT(IN) :: scratch 157 INTEGER, INTENT(IN) :: stdout 158 REAL(DP), INTENT(IN) :: energy_thr, grad_thr, cell_thr 159 LOGICAL, INTENT(IN) :: lmovecell 160 REAL(DP), INTENT(OUT) :: energy_error, grad_error, cell_error 161 LOGICAL, INTENT(OUT) :: step_accepted, stop_bfgs 162 INTEGER, INTENT(OUT) :: istep 163 ! 164 INTEGER :: n, i, j, k, nat 165 LOGICAL :: lwolfe 166 REAL(DP) :: dE0s, den 167 ! ... for scaled coordinates 168 REAL(DP) :: hinv(3,3),g(3,3),ginv(3,3), omega 169 ! 170 ! 171 lwolfe=.false. 172 n = SIZE( pos_in ) + 9 173 nat = size (pos_in) / 3 174 if (nat*3 /= size (pos_in)) call errore('bfgs',' strange dimension',1) 175 ! 176 ! ... work-space allocation 177 ! 178 ALLOCATE( pos( n ) ) 179 ALLOCATE( grad( n ) ) 180 ! 181 ALLOCATE( grad_old( n, bfgs_ndim ) ) 182 ALLOCATE( pos_old( n, bfgs_ndim ) ) 183 ! 184 ALLOCATE( inv_hess( n, n ) ) 185 ! 186 ALLOCATE( pos_p( n ) ) 187 ALLOCATE( grad_p( n ) ) 188 ALLOCATE( step( n ) ) 189 ALLOCATE( step_old( n ) ) 190 ALLOCATE( pos_best( n ) ) 191 ! ... scaled coordinates work-space 192 ALLOCATE( hinv_block( n-9, n-9 ) ) 193 ! ... cell related work-space 194 ALLOCATE( metric( n , n ) ) 195 ! 196 ! ... the BFGS file read (pos & grad) in scaled coordinates 197 ! 198 call invmat(3, h, hinv, omega) 199 ! volume is defined to be positive even for left-handed vector triplet 200 omega = abs(omega) 201 ! 202 hinv_block = 0.d0 203 FORALL ( k=0:nat-1, i=1:3, j=1:3 ) hinv_block(i+3*k,j+3*k) = hinv(i,j) 204 ! 205 ! ... generate metric to work with scaled ionic coordinates 206 g = MATMUL(TRANSPOSE(h),h) 207 call invmat(3,g,ginv) 208 metric = 0.d0 209 FORALL ( k=0:nat-1, i=1:3, j=1:3 ) metric(i+3*k,j+3*k) = g(i,j) 210 FORALL ( k=nat:nat+2, i=1:3, j=1:3 ) metric(i+3*k,j+3*k) = 0.04 * omega * ginv(i,j) 211 ! 212 ! ... generate bfgs vectors for the degrees of freedom and their gradients 213 pos = 0.0 214 pos(1:n-9) = pos_in 215 if (lmovecell) FORALL( i=1:3, j=1:3) pos( n-9 + j+3*(i-1) ) = h(i,j) 216 grad = 0.0 217 grad(1:n-9) = grad_in 218 if (lmovecell) FORALL( i=1:3, j=1:3) grad( n-9 + j+3*(i-1) ) = fcell(i,j)*iforceh(i,j) 219 ! 220 ! if the cell moves the quantity to be minimized is the enthalpy 221 IF ( lmovecell ) fname="enthalpy" 222 ! 223 CALL read_bfgs_file( pos, grad, fixion, energy, scratch, n, stdout ) 224 ! 225 scf_iter = scf_iter + 1 226 istep = scf_iter 227 ! 228 ! ... convergence is checked here 229 ! 230 energy_error = ABS( energy_p - energy ) 231 ! 232 ! ... obscure PGI bug as of v.17.4 233 ! 234#if defined (__PGI) 235 grad_in = MATMUL( TRANSPOSE(hinv_block), grad(1:n-9) ) 236 grad_error = MAXVAL( ABS( grad_in ) ) 237#else 238 grad_error = MAXVAL( ABS( MATMUL( TRANSPOSE(hinv_block), grad(1:n-9)) ) ) 239#endif 240 conv_bfgs = energy_error < energy_thr 241 conv_bfgs = conv_bfgs .AND. ( grad_error < grad_thr ) 242 ! 243 IF( lmovecell) THEN 244 cell_error = MAXVAL( ABS( MATMUL ( TRANSPOSE ( RESHAPE( grad(n-8:n), (/ 3, 3 /) ) ),& 245 TRANSPOSE(h) ) ) ) / omega 246 conv_bfgs = conv_bfgs .AND. ( cell_error < cell_thr ) 247#undef DEBUG 248#if defined(DEBUG) 249 write (*,'(3f15.10)') TRANSPOSE ( RESHAPE( grad(n-8:n), (/ 3, 3 /) ) ) 250 write (*,*) 251 write (*,'(3f15.10)') TRANSPOSE(h) 252 write (*,*) 253 write (*,'(3f15.10)') MATMUL (TRANSPOSE( RESHAPE( grad(n-8:n), (/ 3, 3 /) ) ),& 254 TRANSPOSE(h) ) / omega 255 write (*,*) 256 write (*,*) cell_error/cell_thr*0.5d0 257#endif 258 END IF 259 ! 260 ! ... converged (or useless to go on): quick return 261 ! 262 IF ( .NOT. conv_bfgs .AND. ( tr_min_hit > 1 ) ) CALL infomsg( 'bfgs',& 263 'history already reset at previous step: stopping' ) 264 conv_bfgs = conv_bfgs .OR. ( tr_min_hit > 1 ) 265 ! 266 WRITE(stdout, '(5X,"Energy error",T30,"= ",1PE12.1)') energy_error 267 WRITE(stdout, '(5X,"Gradient error",T30,"= ",1PE12.1)') grad_error 268 IF( lmovecell ) WRITE(stdout, & 269 '(5X,"Cell gradient error",T30,"= ",1PE12.1,/)') cell_error 270 ! 271 IF ( conv_bfgs ) GOTO 1000 272 ! 273 ! ... some output is written 274 ! 275 WRITE( UNIT = stdout, & 276 & FMT = '(/,5X,"number of scf cycles",T30,"= ",I3)' ) scf_iter 277 WRITE( UNIT = stdout, & 278 & FMT = '(5X,"number of bfgs steps",T30,"= ",I3,/)' ) bfgs_iter 279 IF ( scf_iter > 1 ) WRITE( UNIT = stdout, & 280 & FMT = '(5X,A," old",T30,"= ",F18.10," Ry")' ) fname,energy_p 281 WRITE( UNIT = stdout, & 282 & FMT = '(5X,A," new",T30,"= ",F18.10," Ry",/)' ) fname,energy 283 ! 284 ! ... the bfgs algorithm starts here 285 ! 286 IF ( .NOT. energy_wolfe_condition( energy ) .AND. (scf_iter > 1) ) THEN 287 ! 288 ! ... the previous step is rejected, line search goes on 289 ! 290 step_accepted = .FALSE. 291 ! 292 WRITE( UNIT = stdout, & 293 & FMT = '(5X,"CASE: ",A,"_new > ",A,"_old",/)' ) fname,fname 294 ! 295 ! ... the new trust radius is obtained by quadratic interpolation 296 ! 297 ! ... E(s) = a*s*s + b*s + c ( we use E(0), dE(0), E(s') ) 298 ! 299 ! ... s_min = - 0.5*( dE(0)*s'*s' ) / ( E(s') - E(0) - dE(0)*s' ) 300 ! 301 if (abs(scnorm(step_old(:))-1._DP) > 1.d-10) call errore('bfgs', & 302 ' step_old is NOT normalized ',1) 303 ! (normalized) search direction is the same as in previous step 304 step(:) = step_old(:) 305 ! 306 dE0s = ( grad_p(:) .dot. step(:) ) * trust_radius_old 307 IF (dE0s > 0._DP ) CALL errore( 'bfgs', & 308 'dE0s is positive which should never happen', 1 ) 309 den = energy - energy_p - dE0s 310 ! 311 ! estimate new trust radius by interpolation 312 trust_radius = - 0.5_DP*dE0s*trust_radius_old / den 313 ! 314 WRITE( UNIT = stdout, & 315 & FMT = '(5X,"new trust radius",T30,"= ",F18.10," bohr")' ) & 316 trust_radius 317 ! 318 ! ... values from the last successful bfgs step are restored 319 ! 320 pos(:) = pos_p(:) 321 energy = energy_p 322 grad(:) = grad_p(:) 323 ! 324 IF ( trust_radius < trust_radius_min ) THEN 325 ! 326 ! ... the history is reset ( the history can be reset at most two 327 ! ... consecutive times ) 328 ! 329 WRITE( UNIT = stdout, & 330 FMT = '(/,5X,"trust_radius < trust_radius_min")' ) 331 WRITE( UNIT = stdout, FMT = '(/,5X,"resetting bfgs history",/)' ) 332 ! 333 ! ... if tr_min_hit=1 the history has already been reset at the 334 ! ... previous step : something is going wrong 335 ! 336 IF ( tr_min_hit == 1 ) THEN 337 CALL infomsg( 'bfgs', & 338 'history already reset at previous step: stopping' ) 339 tr_min_hit = 2 340 ELSE 341 tr_min_hit = 1 342 END IF 343 ! 344 CALL reset_bfgs( n ) 345 ! 346 step(:) = - ( inv_hess(:,:) .times. grad(:) ) 347 if (lmovecell) FORALL( i=1:3, j=1:3) step( n-9 + j+3*(i-1) ) = step( n-9 + j+3*(i-1) )*iforceh(i,j) 348 ! normalize step but remember its length 349 nr_step_length = scnorm(step) 350 step(:) = step(:) / nr_step_length 351 ! 352 trust_radius = min(trust_radius_ini, nr_step_length) 353 ! 354 ELSE 355 ! 356 tr_min_hit = 0 357 ! 358 END IF 359 ! 360 ELSE 361 ! 362 ! ... a new bfgs step is done 363 ! 364 bfgs_iter = bfgs_iter + 1 365 ! 366 IF ( bfgs_iter == 1 ) THEN 367 ! 368 ! ... first iteration 369 ! 370 step_accepted = .FALSE. 371 ! 372 ELSE 373 ! 374 step_accepted = .TRUE. 375 ! 376 nr_step_length_old = nr_step_length 377 ! 378 WRITE( UNIT = stdout, & 379 & FMT = '(5X,"CASE: ",A,"_new < ",A,"_old",/)' ) fname,fname 380 ! 381 CALL check_wolfe_conditions( lwolfe, energy, grad ) 382 ! 383 CALL update_inverse_hessian( pos, grad, n, stdout ) 384 ! 385 END IF 386 ! compute new search direction and store NR step length 387 IF ( bfgs_ndim > 1 ) THEN 388 ! 389 ! ... GDIIS extrapolation 390 ! 391 CALL gdiis_step() 392 ! 393 ELSE 394 ! 395 ! ... standard Newton-Raphson step 396 ! 397 step(:) = - ( inv_hess(:,:) .times. grad(:) ) 398 if (lmovecell) FORALL( i=1:3, j=1:3) step( n-9 + j+3*(i-1) ) = step( n-9 + j+3*(i-1) )*iforceh(i,j) 399 ! 400 END IF 401 IF ( ( grad(:) .dot. step(:) ) > 0.0_DP ) THEN 402 ! 403 WRITE( UNIT = stdout, & 404 FMT = '(5X,"uphill step: resetting bfgs history",/)' ) 405 ! 406 CALL reset_bfgs( n ) 407 step(:) = - ( inv_hess(:,:) .times. grad(:) ) 408 if (lmovecell) FORALL( i=1:3, j=1:3) step( n-9 + j+3*(i-1) ) = step( n-9 + j+3*(i-1) )*iforceh(i,j) 409 ! 410 END IF 411 ! 412 ! normalize the step and save the step length 413 nr_step_length = scnorm(step) 414 step(:) = step(:) / nr_step_length 415 ! 416 ! ... the new trust radius is computed 417 ! 418 IF ( bfgs_iter == 1 ) THEN 419 ! 420 trust_radius = min(trust_radius_ini, nr_step_length) 421 tr_min_hit = 0 422 ! 423 ELSE 424 ! 425 CALL compute_trust_radius( lwolfe, energy, grad, n, stdout ) 426 ! 427 END IF 428 ! 429 WRITE( UNIT = stdout, & 430 & FMT = '(5X,"new trust radius",T30,"= ",F18.10," bohr")' ) & 431 trust_radius 432 ! 433 END IF 434 ! 435 ! ... step along the bfgs direction 436 ! 437 IF ( nr_step_length < eps16 ) & 438 CALL errore( 'bfgs', 'NR step-length unreasonably short', 1 ) 439 ! 440 ! ... information required by next iteration is saved here ( this must 441 ! ... be done before positions are updated ) 442 ! 443 CALL write_bfgs_file( pos, energy, grad, scratch ) 444 ! 445 ! ... positions and cell are updated 446 ! 447 pos(:) = pos(:) + trust_radius * step(:) 448 ! 4491000 stop_bfgs = conv_bfgs 450 ! ... input ions+cell variables 451 IF ( lmovecell ) FORALL( i=1:3, j=1:3) h(i,j) = pos( n-9 + j+3*(i-1) ) 452 pos_in = pos(1:n-9) 453 ! ... update forces 454 grad_in = grad(1:n-9) 455 ! 456 ! ... work-space deallocation 457 ! 458 DEALLOCATE( pos ) 459 DEALLOCATE( grad ) 460 DEALLOCATE( pos_p ) 461 DEALLOCATE( grad_p ) 462 DEALLOCATE( pos_old ) 463 DEALLOCATE( grad_old ) 464 DEALLOCATE( inv_hess ) 465 DEALLOCATE( step ) 466 DEALLOCATE( step_old ) 467 DEALLOCATE( pos_best ) 468 DEALLOCATE( hinv_block ) 469 DEALLOCATE( metric ) 470 ! 471 RETURN 472 ! 473 CONTAINS 474 ! 475 !-------------------------------------------------------------------- 476 SUBROUTINE gdiis_step() 477 !-------------------------------------------------------------------- 478 USE basic_algebra_routines 479 IMPLICIT NONE 480 ! 481 REAL(DP), ALLOCATABLE :: res(:,:), overlap(:,:), work(:) 482 INTEGER, ALLOCATABLE :: iwork(:) 483 INTEGER :: k, k_m, info 484 REAL(DP) :: gamma0 485 ! 486 ! 487 gdiis_iter = gdiis_iter + 1 488 ! 489 k = MIN( gdiis_iter, bfgs_ndim ) 490 k_m = k + 1 491 ! 492 ALLOCATE( res( n, k ) ) 493 ALLOCATE( overlap( k_m, k_m ) ) 494 ALLOCATE( work( k_m ), iwork( k_m ) ) 495 ! 496 work(:) = 0.0_DP 497 iwork(:) = 0 498 ! 499 ! ... the new direction is added to the workspace 500 ! 501 DO i = bfgs_ndim, 2, -1 502 ! 503 pos_old(:,i) = pos_old(:,i-1) 504 grad_old(:,i) = grad_old(:,i-1) 505 ! 506 END DO 507 ! 508 pos_old(:,1) = pos(:) 509 grad_old(:,1) = grad(:) 510 ! 511 ! ... |res_i> = H^-1 \times |g_i> 512 ! 513 CALL DGEMM( 'N', 'N', n, k, n, 1.0_DP, & 514 inv_hess, n, grad_old, n, 0.0_DP, res, n ) 515 ! 516 ! ... overlap_ij = <grad_i|res_j> 517 ! 518 CALL DGEMM( 'T', 'N', k, k, n, 1.0_DP, & 519 res, n, res, n, 0.0_DP, overlap, k_m ) 520 ! 521 overlap( :, k_m) = 1.0_DP 522 overlap(k_m, : ) = 1.0_DP 523 overlap(k_m,k_m) = 0.0_DP 524 ! 525 ! make sure the overlap matrix is not singular 526 ! 527 FORALL( i = 1:k ) overlap(i,i) = overlap(i,i) * ( 1.0_DP + eps8 ) ; overlap(k_m,k_m) = eps8 528 ! 529 ! ... overlap is inverted via Bunch-Kaufman diagonal pivoting method 530 ! 531 CALL DSYTRF( 'U', k_m, overlap, k_m, iwork, work, k_m, info ) 532 CALL DSYTRI( 'U', k_m, overlap, k_m, iwork, work, info ) 533 CALL errore( 'gdiis_step', 'error in Bunch-Kaufman inversion', info ) 534 ! 535 ! ... overlap is symmetrised 536 ! 537 FORALL( i = 1:k_m, j = 1:k_m, j > i ) overlap(j,i) = overlap(i,j) 538 ! 539 pos_best(:) = 0.0_DP 540 step(:) = 0.0_DP 541 ! 542 DO i = 1, k 543 ! 544 gamma0 = overlap(k_m,i) 545 ! 546 pos_best(:) = pos_best(:) + gamma0*pos_old(:,i) 547 ! 548 step(:) = step(:) - gamma0*res(:,i) 549 ! 550 END DO 551 ! 552 ! ... the step must be consistent with the last positions 553 ! 554 step(:) = step(:) + ( pos_best(:) - pos(:) ) 555 ! 556 IF ( ( grad(:) .dot. step(:) ) > 0.0_DP ) THEN 557 ! 558 ! ... if the extrapolated direction is uphill use only the 559 ! ... last gradient and reset gdiis history 560 ! 561 step(:) = - ( inv_hess(:,:) .times. grad(:) ) 562 if (lmovecell) FORALL( i=1:3, j=1:3) step( n-9 + j+3*(i-1) ) = step( n-9 + j+3*(i-1) )*iforceh(i,j) 563 ! 564 gdiis_iter = 0 565 ! 566 END IF 567 ! 568 DEALLOCATE( res, overlap, work, iwork ) 569 ! 570 END SUBROUTINE gdiis_step 571 ! 572 END SUBROUTINE bfgs 573 ! 574 !------------------------------------------------------------------------ 575 SUBROUTINE reset_bfgs( n ) 576 !------------------------------------------------------------------------ 577 ! ... inv_hess in re-initialized to the initial guess 578 ! ... defined as the inverse metric 579 ! 580 INTEGER, INTENT(IN) :: n 581 ! 582 call invmat(n, metric, inv_hess) 583 ! 584 gdiis_iter = 0 585 ! 586 END SUBROUTINE reset_bfgs 587 ! 588 !------------------------------------------------------------------------ 589 SUBROUTINE read_bfgs_file( pos, grad, fixion, energy, scratch, n, stdout ) 590 !------------------------------------------------------------------------ 591 ! 592 IMPLICIT NONE 593 ! 594 REAL(DP), INTENT(INOUT) :: pos(:) 595 REAL(DP), INTENT(INOUT) :: grad(:) 596 INTEGER, INTENT(IN) :: fixion(:) 597 CHARACTER(LEN=*), INTENT(IN) :: scratch 598 INTEGER, INTENT(IN) :: n 599 INTEGER, INTENT(IN) :: stdout 600 REAL(DP), INTENT(INOUT) :: energy 601 ! 602 CHARACTER(LEN=256) :: bfgs_file 603 LOGICAL :: file_exists 604 ! 605 ! 606 bfgs_file = TRIM( scratch ) // TRIM( prefix ) // '.bfgs' 607 ! 608 INQUIRE( FILE = TRIM( bfgs_file ) , EXIST = file_exists ) 609 ! 610 IF ( file_exists ) THEN 611 ! 612 ! ... bfgs is restarted from file 613 ! 614 OPEN( UNIT = iunbfgs, FILE = TRIM( bfgs_file ), & 615 STATUS = 'UNKNOWN', ACTION = 'READ' ) 616 ! 617 READ( iunbfgs, * ) pos_p 618 READ( iunbfgs, * ) grad_p 619 READ( iunbfgs, * ) scf_iter 620 READ( iunbfgs, * ) bfgs_iter 621 READ( iunbfgs, * ) gdiis_iter 622 READ( iunbfgs, * ) energy_p 623 READ( iunbfgs, * ) pos_old 624 READ( iunbfgs, * ) grad_old 625 READ( iunbfgs, * ) inv_hess 626 READ( iunbfgs, * ) tr_min_hit 627 READ( iunbfgs, * ) nr_step_length 628 ! 629 CLOSE( UNIT = iunbfgs ) 630 ! 631 step_old = ( pos(:) - pos_p(:) ) 632 trust_radius_old = scnorm( step_old ) 633 step_old = step_old / trust_radius_old 634 ! 635 ELSE 636 ! 637 ! ... bfgs initialization 638 ! 639 WRITE( UNIT = stdout, FMT = '(/,5X,"BFGS Geometry Optimization")' ) 640 ! 641 ! initialize the inv_hess to the inverse of the metric 642 call invmat(n, metric, inv_hess) 643 ! 644 pos_p = 0.0_DP 645 grad_p = 0.0_DP 646 scf_iter = 0 647 bfgs_iter = 0 648 gdiis_iter = 0 649 energy_p = energy 650 step_old = 0.0_DP 651 nr_step_length = 0.0_DP 652 ! 653 trust_radius_old = trust_radius_ini 654 ! 655 pos_old(:,:) = 0.0_DP 656 grad_old(:,:) = 0.0_DP 657 ! 658 tr_min_hit = 0 659 ! 660 END IF 661 ! 662 END SUBROUTINE read_bfgs_file 663 ! 664 !------------------------------------------------------------------------ 665 SUBROUTINE write_bfgs_file( pos, energy, grad, scratch ) 666 !------------------------------------------------------------------------ 667 ! 668 IMPLICIT NONE 669 ! 670 REAL(DP), INTENT(IN) :: pos(:) 671 REAL(DP), INTENT(IN) :: energy 672 REAL(DP), INTENT(IN) :: grad(:) 673 CHARACTER(LEN=*), INTENT(IN) :: scratch 674 ! 675 ! 676 OPEN( UNIT = iunbfgs, FILE = TRIM( scratch )//TRIM( prefix )//'.bfgs', & 677 STATUS = 'UNKNOWN', ACTION = 'WRITE' ) 678 ! 679 WRITE( iunbfgs, * ) pos 680 WRITE( iunbfgs, * ) grad 681 WRITE( iunbfgs, * ) scf_iter 682 WRITE( iunbfgs, * ) bfgs_iter 683 WRITE( iunbfgs, * ) gdiis_iter 684 WRITE( iunbfgs, * ) energy 685 WRITE( iunbfgs, * ) pos_old 686 WRITE( iunbfgs, * ) grad_old 687 WRITE( iunbfgs, * ) inv_hess 688 WRITE( iunbfgs, * ) tr_min_hit 689 WRITE( iunbfgs, * ) nr_step_length 690 ! 691 CLOSE( UNIT = iunbfgs ) 692 ! 693 END SUBROUTINE write_bfgs_file 694 ! 695 !------------------------------------------------------------------------ 696 SUBROUTINE update_inverse_hessian( pos, grad, n, stdout ) 697 !------------------------------------------------------------------------ 698 ! 699 IMPLICIT NONE 700 ! 701 REAL(DP), INTENT(IN) :: pos(:) 702 REAL(DP), INTENT(IN) :: grad(:) 703 INTEGER, INTENT(IN) :: n 704 INTEGER, INTENT(IN) :: stdout 705 INTEGER :: info 706 ! 707 REAL(DP), ALLOCATABLE :: y(:), s(:) 708 REAL(DP), ALLOCATABLE :: Hy(:), yH(:) 709 REAL(DP) :: sdoty, sBs, Theta 710 REAL(DP), ALLOCATABLE :: B(:,:) 711 ! 712 ALLOCATE( y( n ), s( n ), Hy( n ), yH( n ) ) 713 ! 714 s(:) = pos(:) - pos_p(:) 715 y(:) = grad(:) - grad_p(:) 716 ! 717 sdoty = ( s(:) .dot. y(:) ) 718 ! 719 IF ( ABS( sdoty ) < eps16 ) THEN 720 ! 721 ! ... the history is reset 722 ! 723 WRITE( stdout, '(/,5X,"WARNING: unexpected ", & 724 & "behaviour in update_inverse_hessian")' ) 725 WRITE( stdout, '( 5X," resetting bfgs history",/)' ) 726 ! 727 CALL reset_bfgs( n ) 728 ! 729 RETURN 730 ! 731 ELSE 732! Conventional Curvature Trap here 733! See section 18.2 (p538-539 ) of Nocedal and Wright "Numerical 734! Optimization"for instance 735! LDM Addition, April 2011 736! 737! While with the Wolfe conditions the Hessian in most cases 738! remains positive definite, if one is far from the minimum 739! and/or "bonds" are being made/broken the curvature condition 740! Hy = s ; or s = By 741! cannot be satisfied if s.y < 0. In addition, if s.y is small 742! compared to s.B.s too greedy a step is taken. 743! 744! The trap below is conventional and "OK", and has been around 745! for ~ 30 years but, unfortunately, is rarely mentioned in 746! introductory texts and hence often neglected. 747! 748! First, solve for inv_hess*t = s ; i.e. t = B*s 749! Use yH as workspace here 750 751 ALLOCATE (B(n,n) ) 752 B = inv_hess 753 yH= s 754 call DPOSV('U',n,1,B,n, yH, n, info) 755! Info .ne. 0 should be trapped ... 756 if(info .ne. 0)write( stdout, '(/,5X,"WARNING: info=",i3," for Hessian")' )info 757 DEALLOCATE ( B ) 758! 759! Calculate s.B.s 760 sBs = ( s(:) .dot. yH(:) ) 761! 762! Now the trap itself 763 if ( sdoty < 0.20D0*sBs ) then 764! Conventional damping 765 Theta = 0.8D0*sBs/(sBs-sdoty) 766 WRITE( stdout, '(/,5X,"WARNING: bfgs curvature condition ", & 767 & "failed, Theta=",F6.3)' )theta 768 y = Theta*y + (1.D0 - Theta)*yH 769 endif 770 END IF 771 ! 772 Hy(:) = ( inv_hess .times. y(:) ) 773 yH(:) = ( y(:) .times. inv_hess ) 774 ! 775 ! ... BFGS update 776 ! 777 inv_hess = inv_hess + 1.0_DP / sdoty * & 778 ( ( 1.0_DP + ( y .dot. Hy ) / sdoty ) * matrix( s, s ) - & 779 ( matrix( s, yH ) + matrix( Hy, s ) ) ) 780 ! 781 DEALLOCATE( y, s, Hy, yH ) 782 ! 783 RETURN 784 ! 785 END SUBROUTINE update_inverse_hessian 786 ! 787 !------------------------------------------------------------------------ 788 SUBROUTINE check_wolfe_conditions( lwolfe, energy, grad ) 789 !------------------------------------------------------------------------ 790 IMPLICIT NONE 791 REAL(DP), INTENT(IN) :: energy 792 REAL(DP), INTENT(IN) :: grad(:) 793 LOGICAL, INTENT(OUT) :: lwolfe 794 ! 795 lwolfe = energy_wolfe_condition ( energy ) .AND. & 796 gradient_wolfe_condition ( grad ) 797 ! 798 END SUBROUTINE check_wolfe_conditions 799 ! 800 !------------------------------------------------------------------------ 801 LOGICAL FUNCTION energy_wolfe_condition ( energy ) 802 !------------------------------------------------------------------------ 803 IMPLICIT NONE 804 REAL(DP), INTENT(IN) :: energy 805 ! 806 energy_wolfe_condition = & 807 ( energy-energy_p ) < w_1 * ( grad_p.dot.step_old ) * trust_radius_old 808 ! 809 END FUNCTION energy_wolfe_condition 810 ! 811 !------------------------------------------------------------------------ 812 LOGICAL FUNCTION gradient_wolfe_condition ( grad ) 813 !------------------------------------------------------------------------ 814 IMPLICIT NONE 815 REAL(DP), INTENT(IN) :: grad(:) 816 ! 817 gradient_wolfe_condition = & 818 ABS( grad .dot. step_old ) < - w_2 * ( grad_p .dot. step_old ) 819 ! 820 END FUNCTION gradient_wolfe_condition 821 ! 822 !------------------------------------------------------------------------ 823 SUBROUTINE compute_trust_radius( lwolfe, energy, grad, n, stdout ) 824 !----------------------------------------------------------------------- 825 ! 826 IMPLICIT NONE 827 ! 828 LOGICAL, INTENT(IN) :: lwolfe 829 REAL(DP), INTENT(IN) :: energy 830 REAL(DP), INTENT(IN) :: grad(:) 831 INTEGER, INTENT(IN) :: n 832 INTEGER, INTENT(IN) :: stdout 833 ! 834 REAL(DP) :: a 835 LOGICAL :: ltest 836 ! 837 ltest = ( energy - energy_p ) < w_1 * ( grad_p .dot. step_old ) * trust_radius_old 838 ! 839 ! The instruction below replaces the original instruction: 840 ! ltest = ltest .AND. ( nr_step_length_old > trust_radius_old ) 841 ! which gives a random result if trust_radius was set equal to 842 ! nr_step_length at previous step. I am not sure what the best 843 ! action should be in that case, though (PG) 844 ! 845 ltest = ltest .AND. ( nr_step_length_old > trust_radius_old + eps8 ) 846 ! 847 IF ( ltest ) THEN 848 a = 1.5_DP 849 ELSE 850 a = 1.1_DP 851 END IF 852 IF ( lwolfe ) a = 2._DP * a 853 ! 854 trust_radius = MIN( trust_radius_max, a*trust_radius_old, nr_step_length ) 855 ! 856 IF ( trust_radius < trust_radius_min ) THEN 857 ! 858 ! ... the history is reset 859 ! 860 ! ... if tr_min_hit the history has already been reset at the 861 ! ... previous step : something is going wrong 862 ! 863 IF ( tr_min_hit == 1 ) THEN 864 CALL infomsg( 'bfgs', & 865 'history already reset at previous step: stopping' ) 866 tr_min_hit = 2 867 ELSE 868 tr_min_hit = 1 869 END IF 870 ! 871 WRITE( UNIT = stdout, & 872 FMT = '(5X,"small trust_radius: resetting bfgs history",/)' ) 873 ! 874 CALL reset_bfgs( n ) 875 step(:) = - ( inv_hess(:,:) .times. grad(:) ) 876 ! 877 nr_step_length = scnorm(step) 878 step(:) = step(:) / nr_step_length 879 ! 880 trust_radius = min(trust_radius_min, nr_step_length ) 881 ! 882 ELSE 883 ! 884 tr_min_hit = 0 885 ! 886 END IF 887 ! 888 END SUBROUTINE compute_trust_radius 889 ! 890 !----------------------------------------------------------------------- 891 REAL(DP) FUNCTION scnorm1( vect ) 892 !----------------------------------------------------------------------- 893 IMPLICIT NONE 894 REAL(DP), INTENT(IN) :: vect(:) 895 ! 896 scnorm1 = SQRT( DOT_PRODUCT( vect , MATMUL( metric, vect ) ) ) 897 ! 898 END FUNCTION scnorm1 899 ! 900 !----------------------------------------------------------------------- 901 REAL(DP) FUNCTION scnorm( vect ) 902 !----------------------------------------------------------------------- 903 IMPLICIT NONE 904 REAL(DP), INTENT(IN) :: vect(:) 905 REAL(DP) :: ss 906 INTEGER :: i,k,l,n 907 ! 908 scnorm = 0._DP 909 n = SIZE (vect) / 3 910 do i=1,n 911 ss = 0._DP 912 do k=1,3 913 do l=1,3 914 ss = ss + & 915 vect(k+(i-1)*3)*metric(k+(i-1)*3,l+(i-1)*3)*vect(l+(i-1)*3) 916 end do 917 end do 918 scnorm = MAX (scnorm, SQRT (ss) ) 919 end do 920 ! 921 END FUNCTION scnorm 922 ! 923 !------------------------------------------------------------------------ 924 SUBROUTINE terminate_bfgs( energy, energy_thr, grad_thr, cell_thr, & 925 lmovecell, stdout, scratch ) 926 !------------------------------------------------------------------------ 927 ! 928 USE io_files, ONLY : delete_if_present 929 ! 930 IMPLICIT NONE 931 REAL(DP), INTENT(IN) :: energy, energy_thr, grad_thr, cell_thr 932 LOGICAL, INTENT(IN) :: lmovecell 933 INTEGER, INTENT(IN) :: stdout 934 CHARACTER(LEN=*), INTENT(IN) :: scratch 935 ! 936 IF ( conv_bfgs ) THEN 937 ! 938 WRITE( UNIT = stdout, & 939 & FMT = '(/,5X,"bfgs converged in ",I3," scf cycles and ", & 940 & I3," bfgs steps")' ) scf_iter, bfgs_iter 941 IF ( lmovecell ) THEN 942 WRITE( UNIT = stdout, & 943 & FMT = '(5X,"(criteria: energy < ",ES8.1," Ry, force < ",ES8.1,& 944 & "Ry/Bohr, cell < ",ES8.1,"kbar)")') energy_thr, grad_thr, cell_thr 945 ELSE 946 WRITE( UNIT = stdout, & 947 & FMT = '(5X,"(criteria: energy < ",ES8.1," Ry, force < ",ES8.1,& 948 & " Ry/Bohr)")') energy_thr, grad_thr 949 END IF 950 WRITE( UNIT = stdout, & 951 & FMT = '(/,5X,"End of BFGS Geometry Optimization")' ) 952 WRITE( UNIT = stdout, & 953 & FMT = '(/,5X,"Final ",A," = ",F18.10," Ry")' ) fname, energy 954 ! 955 CALL delete_if_present( TRIM( scratch ) // TRIM( prefix ) // '.bfgs' ) 956 ! 957 ELSE 958 ! 959 WRITE( UNIT = stdout, & 960 FMT = '(/,5X,"The maximum number of steps has been reached.")' ) 961 WRITE( UNIT = stdout, & 962 FMT = '(/,5X,"End of BFGS Geometry Optimization")' ) 963 ! 964 END IF 965 ! 966 END SUBROUTINE terminate_bfgs 967 ! 968 FUNCTION bfgs_get_n_iter (what) RESULT(n_iter) 969 ! 970 IMPLICIT NONE 971 INTEGER :: n_iter 972 CHARACTER(10),INTENT(IN) :: what 973 SELECT CASE (TRIM(what)) 974 CASE ('bfgs_iter') 975 n_iter = bfgs_iter 976 CASE ( 'scf_iter') 977 n_iter = scf_iter 978 CASE default 979 n_iter = -1 980 END SELECT 981 END FUNCTION bfgs_get_n_iter 982END MODULE bfgs_module 983