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