1 subroutine lbfgd(emin,ggtol,nsd,x,g,diag,w) 2 implicit real (a-h,o-z), integer (i-n) 3 common /lb3/gtol,stpmin,stpmax 4 common /athlp/ iatoms, mxnat 5 common /rdwr/ iun1,iun2,iun3,iun4,iun5 6 logical osingl,dolbfgs,oqscal 7 common /opts/ idebug,imon,iarc,ilog,iout,osingl,dolbfgs,oqscal 8 logical usecut,usesw 9 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 10 logical box,cell,fast 11 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 12 common /cell/ xa,ya,yb,za,zb,zc 13 common /convg/ icyc 14 dimension x(*),g(*),diag(*),w(*) 15 16 write(iun5,*) 'L-BFGS method' 17 write(iun5,*) ' ' 18 19 gtol = 0.9e0 20 stpmin = 1.0e-20 21 stpmax = 1.0e+20 22 iwrite = 1 23 24 n = iatoms*3 25 26 if (cell) then 27 x(n+1) = xa 28 x(n+2) = ya 29 x(n+3) = yb 30 x(n+4) = za 31 x(n+5) = zb 32 x(n+6) = zc 33 34 n = n + 6 35 endif 36 37 m = 5 38 39 40 eps = ggtol 41 42 xtol = 1.0e-16 43 write(iun5,*) 'Gradient tolerance: ',eps 44 write(iun5,*) 'Maximum iterations: ',nsd 45 write(iun5,*) ' ' 46 47 ncyc = 0 48 icyc = 0 49 50 iflag = 0 51 52 do while (.true.) 53 54 call enegrd(emin,x,g) 55 call rmsg(n,g,rmsgrd) 56 57 if (osingl) then 58 write(iun5,'(a,f12.3)') 'Estat=',emin 59 return 60 endif 61 62 do i=1,n 63 g(i) = -g(i) 64 end do 65 66 call dlbfgs(n,m,x,emin,g,diag,eps,xtol,w,iflag) 67 if (iflag.eq.0) return 68 if (iflag.lt.0) then 69 return 70 endif 71 72c We have found a lower energy. 73 74 ncyc = ncyc + 1 75 76c Write out the minimized structure 77 78 if (mod(ncyc,iwrite) .eq. 0) then 79 80c arc file 81 if (iarc.eq.1) call wrtout(iun4,emin) 82 83c intermediate file 84 85 icyc = icyc + 1 86 if (imon.ne.0) call wrmon(icyc,emin) 87 88 if (usesw) then 89 if (mod(icyc,10).eq.0) call bldlst 90 endif 91 92 end if 93 94 if (ncyc.gt.nsd) then 95 write(iun5,*) 'Exceeded maximum interations ',nsd 96 return 97 endif 98 99 end do 100 101 return 102 end 103 104 subroutine dlbfgs(n,m,x,f,g,diag,eps,xtol,w,iflag) 105 implicit real (a-h,o-z), integer (i-n) 106 logical box,cell,fast 107 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 108 dimension x(*),g(*),diag(*),w(*) 109 110c limited memory bfgs method for large scale optimization 111c jorge nocedal 112 113 common /lb3/gtol,stpmin,stpmax 114 common /rdwr/ iun1,iun2,iun3,iun4,iun5 115 116 integer bound,cp,point 117 logical finish 118 save 119 120c initialize 121 122 if (iflag.eq.1) goto 172 123 124 iter = 0 125 126 if (n.le.0.or.m.le.0) goto 196 127 128 if (gtol.le.1.e-04) then 129 write(iun5,245) 130 gtol = 9.e-01 131 endif 132 133 nfun = 1 134 point = 0 135 finish = .false. 136 137 do i=1,n 138 diag(i) = 1.0e0 139 end do 140 141 ispt = n+2*m 142 iypt = ispt+n*m 143 144 do i=1,n 145 w(ispt+i) = -g(i)*diag(i) 146 end do 147 148 gnorm = sqrt(ddot(n,g,g)) 149 gpnorm = sqrt(ddot(n,g,g)/dble(n)) 150 stp1 = 1.0e0/gnorm 151 152c parameters for line search routine 153 154 ftol = 1.0e-4 155 maxfev = 20 156 157 call prtene(nfun,f,gpnorm) 158 159c -------------------- 160c main iteration loop 161c -------------------- 162 163100 continue 164 165 iter = iter+1 166 info = 0 167 bound = iter-1 168 169 if (iter.ne.1) then 170 if (iter.gt.m) bound = m 171 172 ys = ddot(n,w(iypt+npt+1),w(ispt+npt+1)) 173 yy = ddot(n,w(iypt+npt+1),w(iypt+npt+1)) 174 175 do i=1,n 176 diag(i) = ys/yy 177 end do 178 179 cp = point 180 if (point.eq.0) cp = m 181 w(n+cp) = 1.0e0/ys 182 183 do i=1,n 184 w(i) = -g(i) 185 end do 186 187 cp = point 188 189 do i=1,bound 190 cp = cp-1 191 if (cp.eq. -1) cp = m-1 192 sq = ddot(n,w(ispt+cp*n+1),w) 193 inmc = n + m + cp + 1 194 iycn = iypt + cp*n 195 w(inmc) = w(n+cp+1)*sq 196 call daxpy(n,-w(inmc),w(iycn+1),w) 197 end do 198 199 do i=1,n 200 w(i) = diag(i)*w(i) 201 end do 202 203 do i=1,bound 204 205 yr = ddot(n,w(iypt+cp*n+1),w) 206 beta = w(n+cp+1)*yr 207 inmc = n + m + cp + 1 208 beta = w(inmc)-beta 209 iscn = ispt + cp*n 210 211 call daxpy(n,beta,w(iscn+1),w) 212 213 cp = cp + 1 214 if (cp.eq.m) cp = 0 215 216 end do 217 218 do i=1,n 219 w(ispt+point*n+i) = w(i) 220 end do 221 222 endif 223 224 nfev = 0 225 stp = 1.0e0 226 if (iter.eq.1) stp = stp1 227 228 do i=1,n 229 w(i) = g(i) 230 end do 231 232 172 continue 233 234 call mcsrch(n,x,f,g,w(ispt+point*n+1),stp,ftol, 235 & xtol,maxfev,info,nfev,diag) 236 if (info .eq. -1) then 237 iflag = 1 238 return 239 endif 240 241 if (info.ne.1) goto 190 242 243 nfun = nfun + nfev 244 245c compute the new step and gradient change 246 247 npt = point*n 248 do i=1,n 249 w(ispt+npt+i) = stp*w(ispt+npt+i) 250 w(iypt+npt+i) = g(i)-w(i) 251 end do 252 253 point = point + 1 254 if (point.eq.m) point = 0 255 256c termination test 257 258 gnorm = sqrt(ddot(n,g,g)) 259 gpnorm = sqrt(ddot(n,g,g)/dble(n)) 260 xnorm = sqrt(ddot(n,x,x)) 261 xnorm = max(1.0e0,xnorm) 262 263 if (gpnorm .le. eps) finish = .true. 264 265 call prtene(nfun,f,gpnorm) 266 267 if (finish) then 268 write(iun5,*) 'Optimisation Converged' 269 iflag = 0 270 return 271 endif 272 273 goto 100 274 275c end of main iteration loop. error exits. 276 277 190 iflag=-1 278 write(iun5,200) info 279 return 280 281 195 iflag=-2 282 write(iun5,235) i 283 return 284 285 196 iflag= -3 286 write(iun5,240) 287 288c formats 289 290 200 format(/' iflag= -1 ',/' line search failed. see', 291 & ' documentation of routine mcsrch',/' error return', 292 & ' of line search: info= ',i2,/ 293 & ' possible causes: function or gradient are incorrect',/, 294 & ' or incorrect tolerances') 295 235 format(/' iflag= -2',/' the',i5,'-th diagonal element of the',/, 296 & ' inverse hessian approximation is not positive') 297 240 format(/' iflag= -3',/' improper input parameters (n or m', 298 & ' are not positive)') 299 245 format(/' gtol is less than or equal to 1.e-04', 300 & / ' it has been reset to 9.e-01') 301 302 return 303 end 304 305 subroutine daxpy(n,da,dx,dy) 306 implicit real (a-h,o-z), integer (i-n) 307 dimension dx(*),dy(*) 308 309 if (n.le.0.or.da.eq.0.0e0) return 310 311 m = mod(n,4) 312 313 if (m.ne.0) then 314 do i= 1,m 315 dy(i) = dy(i) + da*dx(i) 316 end do 317 318 if (n.lt.4) return 319 320 endif 321 322 mp1 = m + 1 323 324 do i= mp1,n,4 325 dy(i) = dy(i) + da*dx(i) 326 dy(i+1) = dy(i+1) + da*dx(i+1) 327 dy(i+2) = dy(i+2) + da*dx(i+2) 328 dy(i+3) = dy(i+3) + da*dx(i+3) 329 end do 330 331 return 332 end 333 334 real function ddot(n,dx,dy) 335 336c forms the dot product of two vectors. 337 338 implicit real (a-h,o-z), integer (i-n) 339 dimension dx(*),dy(*) 340 341 ddot = 0.0e0 342 dtemp = 0.0e0 343 344 if (n.le.0) return 345 346 347 m = mod(n,5) 348 if (m.ne.0) then 349 350 do i=1,m 351 dtemp = dtemp + dx(i)*dy(i) 352 end do 353 354 if (n.lt.5) goto 60 355 endif 356 357 mp1 = m + 1 358 359 do i = mp1,n,5 360 dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + 361 & dx(i+2)*dy(i+2) + dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4) 362 end do 363 364 36560 ddot = dtemp 366 367 return 368 end 369 370 371 subroutine mcsrch(n,x,f,g,s,stp,ftol,xtol,maxfev,info,nfev,wa) 372 373c line search routine 374 375 implicit real (a-h,o-z), integer (i-n) 376 common /lb3/gtol,stpmin,stpmax 377 common /rdwr/ iun1,iun2,iun3,iun4,iun5 378 save 379 logical brackt,stage1 380 381 dimension x(*),g(*),s(*),wa(*) 382 383 if (info.ne.-1) then 384 385 infoc = 1 386 387 388 if (n .le. 0 .or. stp .le. 0.0e0 .or. ftol .lt. 0.0e0 .or. 389 & gtol .lt. 0.0e0 .or. xtol .lt. 0.0e0 .or. stpmin .lt. 0.0e0 390 & .or. stpmax .lt. stpmin .or. maxfev .le. 0) return 391 392 dginit = 0.0e0 393 394 do j = 1, n 395 dginit = dginit + g(j)*s(j) 396 end do 397 398 if (dginit .ge. 0.0e0) then 399 write(iun5,'(a)') 400 & ' the search direction is not a descent direction' 401 return 402 endif 403 404 405 brackt = .false. 406 stage1 = .true. 407 nfev = 0 408 finit = f 409 dgtest = ftol*dginit 410 width = stpmax - stpmin 411 width1 = width/0.5e0 412 413 do j = 1, n 414 wa(j) = x(j) 415 end do 416 417 stx = 0.0e0 418 fx = finit 419 dgx = dginit 420 sty = 0.0e0 421 fy = finit 422 dgy = dginit 423 424 endif 425 426c start of iteration. 427 428 do while (.true.) 429 430 if (info.ne.-1) then 431 if (brackt) then 432 stmin = min(stx,sty) 433 stmax = max(stx,sty) 434 else 435 stmin = stx 436 stmax = stp + 4.0e0*(stp - stx) 437 end if 438 439 stp = max(stp,stpmin) 440 stp = min(stp,stpmax) 441 442 443 if ((brackt .and. (stp.le.stmin .or. stp.ge.stmax)) 444 & .or. nfev.ge.maxfev-1 .or. infoc.eq.0 445 & .or. (brackt .and. stmax-stmin.le.xtol*stmax)) stp = stx 446 447 do j = 1, n 448 x(j) = wa(j) + stp*s(j) 449 end do 450 451 info = -1 452 return 453 454 endif 455 456 info = 0 457 nfev = nfev + 1 458 dg = 0.0e0 459 460 do j = 1, n 461 dg = dg + g(j)*s(j) 462 end do 463 464 ftest1 = finit + stp*dgtest 465 466c test for convergence. 467 468 if ((brackt .and. (stp.le.stmin .or. stp.ge.stmax)) 469 & .or. infoc.eq.0) info = 6 470 471 if (stp.eq.stpmax .and. 472 & f.le.ftest1 .and. dg.le.dgtest) info = 5 473 474 if (stp.eq.stpmin .and. 475 & (f.gt.ftest1 .or. dg.ge.dgtest)) info = 4 476 477 if (nfev.ge.maxfev) info = 3 478 479 if (brackt .and. stmax-stmin.le.xtol*stmax) info = 2 480 481 if (f.le.ftest1 .and. abs(dg).le.gtol*(-dginit)) info = 1 482 483c check for termination. 484 485 if (info.ne.0) return 486 487 488 if (stage1 .and. f.le.ftest1 .and. 489 & dg.ge.min(ftol,gtol)*dginit) stage1 = .false. 490 491 if (stage1 .and. f.le.fx .and. f.gt.ftest1) then 492 493 fm = f - stp*dgtest 494 fxm = fx - stx*dgtest 495 fym = fy - sty*dgtest 496 dgm = dg - dgtest 497 dgxm = dgx - dgtest 498 dgym = dgy - dgtest 499 500 call mcstep(stx,fxm,dgxm,sty,fym,dgym,stp,fm,dgm, 501 & brackt,stmin,stmax,infoc) 502 503 fx = fxm + stx*dgtest 504 fy = fym + sty*dgtest 505 dgx = dgxm + dgtest 506 dgy = dgym + dgtest 507 508 else 509 510 call mcstep(stx,fx,dgx,sty,fy,dgy,stp,f,dg, 511 & brackt,stmin,stmax,infoc) 512 end if 513 514 if (brackt) then 515 516 if (abs(sty-stx).ge.0.66e0*width1) 517 & stp = stx + 0.5e0*(sty - stx) 518 width1 = width 519 width = abs(sty-stx) 520 521 end if 522 523 end do 524 525 return 526 end 527 528 subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, 529 & stpmin,stpmax,info) 530 implicit real (a-h,o-z), integer (i-n) 531 logical brackt,bound 532 533c the purpose of mcstep is to compute a safeguarded step for 534c a linesearch and to update an interval of uncertainty for 535c a minimizer of the function. 536 537 538 info = 0 539 540c check the input parameters for errors. 541 542 if ((brackt .and. (stp .le. min(stx,sty) .or. 543 & stp .ge. max(stx,sty))) .or. 544 & dx*(stp-stx) .ge. 0.0e0 .or. stpmax .lt. stpmin) return 545 546c determine if the derivatives have opposite sign. 547 548 sgnd = dp*(dx/abs(dx)) 549 550c first case. a higher function value. 551c the minimum is bracketed. if the cubic step is closer 552c to stx than the quadratic step, the cubic step is taken, 553c else the average of the cubic and quadratic steps is taken. 554 555 if (fp .gt. fx) then 556 557 info = 1 558 bound = .true. 559 theta = 3*(fx - fp)/(stp - stx) + dx + dp 560 s = max(abs(theta),abs(dx),abs(dp)) 561 gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s)) 562 if (stp .lt. stx) gamma = -gamma 563 p = (gamma - dx) + theta 564 q = ((gamma - dx) + gamma) + dp 565 r = p/q 566 stpc = stx + r*(stp - stx) 567 stpq = stx + ((dx/((fx-fp)/(stp-stx)+dx))/2)*(stp - stx) 568 569 if (abs(stpc-stx) .lt. abs(stpq-stx)) then 570 stpf = stpc 571 else 572 stpf = stpc + (stpq - stpc)/2 573 end if 574 575 brackt = .true. 576 577c second case. a lower function value and derivatives of 578c opposite sign. the minimum is bracketed. if the cubic 579c step is closer to stx than the quadratic (secant) step, 580c the cubic step is taken, else the quadratic step is taken. 581 582 else if (sgnd .lt. 0.0) then 583 584 info = 2 585 bound = .false. 586 theta = 3*(fx - fp)/(stp - stx) + dx + dp 587 s = max(abs(theta),abs(dx),abs(dp)) 588 gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s)) 589 if (stp .gt. stx) gamma = -gamma 590 p = (gamma - dp) + theta 591 q = ((gamma - dp) + gamma) + dx 592 r = p/q 593 stpc = stp + r*(stx - stp) 594 stpq = stp + (dp/(dp-dx))*(stx - stp) 595 596 if (abs(stpc-stp) .gt. abs(stpq-stp)) then 597 stpf = stpc 598 else 599 stpf = stpq 600 end if 601 602 brackt = .true. 603 604c third case. a lower function value, derivatives of the 605c same sign, and the magnitude of the derivative decreases. 606c the cubic step is only used if the cubic tends to infinity 607c in the direction of the step or if the minimum of the cubic 608c is beyond stp. otherwise the cubic step is defined to be 609c either stpmin or stpmax. the quadratic (secant) step is also 610c computed and if the minimum is bracketed then the the step 611c closest to stx is taken, else the step farthest away is taken. 612 613 else if (abs(dp) .lt. abs(dx)) then 614 615 info = 3 616 bound = .true. 617 theta = 3*(fx - fp)/(stp - stx) + dx + dp 618 s = max(abs(theta),abs(dx),abs(dp)) 619 620c the case gamma = 0 only arises if the cubic does not tend 621c to infinity in the direction of the step. 622 623 gamma = s*sqrt(max(0.0e0,(theta/s)**2 - (dx/s)*(dp/s))) 624 625 if (stp .gt. stx) gamma = -gamma 626 627 p = (gamma - dp) + theta 628 q = (gamma + (dx - dp)) + gamma 629 r = p/q 630 631 if (r .lt. 0.0e0 .and. gamma .ne. 0.0e0) then 632 stpc = stp + r*(stx - stp) 633 else if (stp .gt. stx) then 634 stpc = stpmax 635 else 636 stpc = stpmin 637 endif 638 639 stpq = stp + (dp/(dp-dx))*(stx - stp) 640 641 if (brackt) then 642 643 if (abs(stp-stpc) .lt. abs(stp-stpq)) then 644 stpf = stpc 645 else 646 stpf = stpq 647 end if 648 649 else 650 651 if (abs(stp-stpc) .gt. abs(stp-stpq)) then 652 stpf = stpc 653 else 654 stpf = stpq 655 end if 656 657 end if 658 659c fourth case. a lower function value, derivatives of the 660c same sign, and the magnitude of the derivative does 661c not decrease. if the minimum is not bracketed, the step 662c is either stpmin or stpmax, else the cubic step is taken. 663 664 else 665 666 info = 4 667 bound = .false. 668 669 if (brackt) then 670 theta = 3*(fp - fy)/(sty - stp) + dy + dp 671 s = max(abs(theta),abs(dy),abs(dp)) 672 gamma = s*sqrt((theta/s)**2 - (dy/s)*(dp/s)) 673 if (stp .gt. sty) gamma = -gamma 674 p = (gamma - dp) + theta 675 q = ((gamma - dp) + gamma) + dy 676 r = p/q 677 stpc = stp + r*(sty - stp) 678 stpf = stpc 679 else if (stp .gt. stx) then 680 stpf = stpmax 681 else 682 stpf = stpmin 683 end if 684 685 endif 686 687c update the interval of uncertainty. this update does not 688c depend on the new step or the case analysis above. 689 690 if (fp .gt. fx) then 691 692 sty = stp 693 fy = fp 694 dy = dp 695 696 else 697 698 if (sgnd.lt.0.0e0) then 699 sty = stx 700 fy = fx 701 dy = dx 702 end if 703 704 stx = stp 705 fx = fp 706 dx = dp 707 708 endif 709 710c compute the new step and safeguard it. 711 712 stpf = min(stpmax,stpf) 713 stpf = max(stpmin,stpf) 714 stp = stpf 715 716 if (brackt .and. bound) then 717 718 stt = stx+0.66*(sty-stx) 719 720 if (sty .gt. stx) then 721 stp = min(stt,stp) 722 else 723 stp = max(stt,stp) 724 endif 725 726 endif 727 728 return 729 end 730 731 subroutine rmsg(n,g1,rmsgrd) 732 implicit real (a-h,o-z), integer (i-n) 733 dimension g1(*) 734 735 g1g1 = 0.0e0 736 do i=1,n 737 g1g1 = g1g1 + g1(i)*g1(i) 738 end do 739 740 rmsgrd = sqrt(g1g1/dble(n)) 741 742 return 743 end 744 745