+
+!
+!
+! ###################################################
+! ## COPYRIGHT (C) 1999 by Jay William Ponder ##
+! ## All Rights Reserved ##
+! ###################################################
+!
+! ##############################################################
+! ## ##
+! ## subroutine lbfgs -- limited memory BFGS optimization ##
+! ## ##
+! ##############################################################
+!
+!
+! "lbfgs" is a limited memory BFGS quasi-newton nonlinear
+! optimization routine
+!
+! literature references:
+!
+! J. Nocedal, "Updating Quasi-Newton Matrices with Limited
+! Storage", Mathematics of Computation, 35, 773-782 (1980)
+!
+! D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
+! for Large Scale Optimization", Mathematical Programming,
+! 45, 503-528 (1989)
+!
+! J. Nocedal and S. J. Wright, "Numerical Optimization",
+! Springer-Verlag, New York, 1999, Section 9.1
+!
+! variables and parameters:
+!
+! nvar number of parameters in the objective function
+! x0 contains starting point upon input, upon return
+! contains the best point found
+! minimum during optimization contains best current function
+! value; returns final best function value
+! grdmin normal exit if rms gradient gets below this value
+! ncalls total number of function/gradient evaluations
+!
+! required external routines:
+!
+! fgvalue function to evaluate function and gradient values
+! optsave subroutine to write out info about current status
+!
+!
+ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
+ use inform
+ use iounit
+ use keys
+ use linmin
+ use math2
+ use minima
+ use output
+ use scales
+ use MD_data, only: itime_mat
+#ifdef LBFGS
+ use energy_data,only:statusbf,niter,ncalls
+#endif
+ implicit none
+ integer i,j,k,m
+ integer nvar,next
+ integer msav,muse
+#ifndef LBFGS
+ integer niter,ncalls
+ character*9 statusbf
+#endif
+ integer nerr,maxerr
+ real*8 f,f_old,fgvalue
+ real*8 f_move,x_move
+ real*8 g_norm,g_rms
+ real*8 minimum,grdmin
+ real*8 angle,rms,beta
+ real*8 ys,yy,gamma
+ real*8 x0(*)
+ real*8, allocatable :: rho(:)
+ real*8, allocatable :: alpha(:)
+ real*8, allocatable :: x_old(:)
+ real*8, allocatable :: g(:)
+ real*8, allocatable :: g_old(:)
+ real*8, allocatable :: p(:)
+ real*8, allocatable :: q(:)
+ real*8, allocatable :: r(:)
+ real*8, allocatable :: h0(:)
+ real*8, allocatable :: s(:,:)
+ real*8, allocatable :: y(:,:)
+ logical done
+ character*9 blank
+ character*20 keyword
+ character*240 record
+ character*240 string
+ external fgvalue,optsave
+! common /lbfgstat/ status,niter,ncalls
+!c
+!c
+!c initialize some values to be used below
+!c
+ ncalls = 0
+ rms = sqrt(dble(nvar))
+ if (coordtype .eq. 'CARTESIAN') then
+ rms = rms / sqrt(3.0d0)
+ else if (coordtype .eq. 'RIGIDBODY') then
+ rms = rms / sqrt(6.0d0)
+ end if
+ blank = ' '
+ done = .false.
+ nerr = 0
+ maxerr = 2
+!c
+!c perform dynamic allocation of some global arrays
+!c
+ if (.not. allocated(scale)) allocate (scale(nvar))
+!c
+!c set default values for variable scale factors
+!c
+ if (.not. set_scale) then
+ do i = 1, nvar
+ if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0
+ end do
+ end if
+!c
+!c set default parameters for the optimization
+!c
+ msav = min(nvar,20)
+ if (fctmin .eq. 0.0d0) fctmin = -100000000.0d0
+ if (maxiter .eq. 0) maxiter = 1000000
+ if (nextiter .eq. 0) nextiter = 1
+ if (jprint .lt. 0) jprint = 1
+ if (iwrite .lt. 0) iwrite = 1
+ write(iout,*) "maxiter",maxiter
+!c
+!c set default parameters for the line search
+!c
+ if (stpmax .eq. 0.0d0) stpmax = 5.0d0
+ stpmin = 1.0d-16
+ cappa = 0.9d0
+ slpmax = 10000.0d0
+ angmax = 180.0d0
+ intmax = 5
+!c
+!c search the keywords for optimization parameters
+!c
+#ifdef LBFGSREAD
+ do i = 1, nkey
+ next = 1
+ record = keyline(i)
+ call gettext (record,keyword,next)
+ call upcase (keyword)
+ string = record(next:240)
+ if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
+ read (string,*,err=10,end=10) msav
+ msav = max(0,min(msav,nvar))
+ else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
+ msav = 0
+ else if (keyword(1:7) .eq. 'FCTMIN ') then
+ read (string,*,err=10,end=10) fctmin
+ else if (keyword(1:8) .eq. 'MAXITER ') then
+ read (string,*,err=10,end=10) maxiter
+ else if (keyword(1:8) .eq. 'STEPMAX ') then
+ read (string,*,err=10,end=10) stpmax
+ else if (keyword(1:8) .eq. 'STEPMIN ') then
+ read (string,*,err=10,end=10) stpmin
+ else if (keyword(1:6) .eq. 'CAPPA ') then
+ read (string,*,err=10,end=10) cappa
+ else if (keyword(1:9) .eq. 'SLOPEMAX ') then
+ read (string,*,err=10,end=10) slpmax
+ else if (keyword(1:7) .eq. 'ANGMAX ') then
+ read (string,*,err=10,end=10) angmax
+ else if (keyword(1:7) .eq. 'INTMAX ') then
+ read (string,*,err=10,end=10) intmax
+ end if
+ 10 continue
+ end do
+#endif
+!c
+!c print header information about the optimization method
+!c
+ if (jprint .gt. 0) then
+ if (msav .eq. 0) then
+ write (jout,20)
+ 20 format (/,' Steepest Descent Gradient Optimization :')
+ write (jout,30)
+ 30 format (/,' SD Iter F Value G RMS F Move',&
+ ' X Move Angle FG Call Comment',/)
+ else
+ write (jout,40)
+ 40 format (/,' Limited Memory BFGS Quasi-Newton',&
+ ' Optimization :')
+ write (jout,50)
+ 50 format (/,' QN Iter F Value G RMS F Move',&
+ ' X Move Angle FG Call Comment',/)
+ end if
+ flush (jout)
+ end if
+!c
+!c perform dynamic allocation of some local arrays
+!c
+ allocate (x_old(nvar))
+ allocate (g(nvar))
+ allocate (g_old(nvar))
+ allocate (p(nvar))
+ allocate (q(nvar))
+ allocate (r(nvar))
+ allocate (h0(nvar))
+ if (msav .ne. 0) then
+ allocate (rho(msav))
+ allocate (alpha(msav))
+ allocate (s(nvar,msav))
+ allocate (y(nvar,msav))
+ end if
+!c
+!c evaluate the function and get the initial gradient
+!c
+ niter = nextiter - 1
+ maxiter = niter + maxiter
+ ncalls = ncalls + 1
+! itime_mat=niter
+ f = fgvalue (x0,g)
+ f_old = f
+ m = 0
+ gamma = 1.0d0
+ g_norm = 0.0d0
+ g_rms = 0.0d0
+ do i = 1, nvar
+ g_norm = g_norm + g(i)*g(i)
+ g_rms = g_rms + (g(i)*scale(i))**2
+ end do
+ g_norm = sqrt(g_norm)
+ g_rms = sqrt(g_rms) / rms
+ f_move = 0.5d0 * stpmax * g_norm
+!c
+!c print initial information prior to first iteration
+!c
+ write(jout,*) "before first"
+ if (jprint .gt. 0) then
+ if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
+ write (jout,60) niter,f,g_rms,ncalls
+ 60 format (i6,f14.4,f11.4,2x,i7)
+ else
+ write (jout,70) niter,f,g_rms,ncalls
+ 70 format (i6,d14.4,d11.4,2x,i7)
+ end if
+ flush (jout)
+ end if
+!c
+!c write initial intermediate prior to first iteration
+!c
+ if (iwrite .gt. 0) call optsave (niter,f,x0)
+!c
+!c tests of the various termination criteria
+!c
+ if (niter .ge. maxiter) then
+ statusbf = 'IterLimit'
+ done = .true.
+ end if
+ if (f .le. fctmin) then
+ statusbf = 'SmallFct '
+ done = .true.
+ end if
+ if (g_rms .le. grdmin) then
+ statusbf = 'SmallGrad'
+ done = .true.
+ end if
+!c
+!c start of a new limited memory BFGS iteration
+!c
+ do while (.not. done)
+ niter = niter + 1
+!c write (jout,*) "LBFGS niter",niter
+ muse = min(niter-1,msav)
+ m = m + 1
+ if (m .gt. msav) m = 1
+!c
+!c estimate Hessian diagonal and compute the Hg product
+!c
+ do i = 1, nvar
+ h0(i) = gamma
+ q(i) = g(i)
+ end do
+ k = m
+ do j = 1, muse
+ k = k - 1
+ if (k .eq. 0) k = msav
+ alpha(k) = 0.0d0
+ do i = 1, nvar
+ alpha(k) = alpha(k) + s(i,k)*q(i)
+ end do
+ alpha(k) = alpha(k) * rho(k)
+ do i = 1, nvar
+ q(i) = q(i) - alpha(k)*y(i,k)
+ end do
+ end do
+ do i = 1, nvar
+ r(i) = h0(i) * q(i)
+ end do
+ do j = 1, muse
+ beta = 0.0d0
+ do i = 1, nvar
+ beta = beta + y(i,k)*r(i)
+ end do
+ beta = beta * rho(k)
+ do i = 1, nvar
+ r(i) = r(i) + s(i,k)*(alpha(k)-beta)
+ end do
+ k = k + 1
+ if (k .gt. msav) k = 1
+ end do
+!c
+!c set search direction and store current point and gradient
+!c
+ do i = 1, nvar
+ p(i) = -r(i)
+ x_old(i) = x0(i)
+ g_old(i) = g(i)
+ end do
+!c
+!c perform line search along the new conjugate direction
+!c
+ statusbf = blank
+!c write (jout,*) "Before search"
+ call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,statusbf)
+!c write (jout,*) "After search"
+!c
+!c update variables based on results of this iteration
+!c
+ if (msav .ne. 0) then
+ ys = 0.0d0
+ yy = 0.0d0
+ do i = 1, nvar
+ s(i,m) = x0(i) - x_old(i)
+ y(i,m) = g(i) - g_old(i)
+ ys = ys + y(i,m)*s(i,m)
+ yy = yy + y(i,m)*y(i,m)
+ end do
+ gamma = abs(ys/yy)
+ rho(m) = 1.0d0 / ys
+ end if
+!c
+!c get the sizes of the moves made during this iteration
+!c
+ f_move = f_old - f
+! if (f_move.eq.0) f_move=1.0d0
+ f_old = f
+ x_move = 0.0d0
+ do i = 1, nvar
+ x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
+ end do
+ x_move = sqrt(x_move) / rms
+ if (coordtype .eq. 'INTERNAL') then
+ x_move = radian * x_move
+ end if
+!c
+!c compute the rms gradient per optimization parameter
+!c
+ g_rms = 0.0d0
+ do i = 1, nvar
+ g_rms = g_rms + (g(i)*scale(i))**2
+ end do
+ g_rms = sqrt(g_rms) / rms
+!c
+!c test for error due to line search problems
+!c
+ if (statusbf.eq.'BadIntpln' .or. statusbf.eq.'IntplnErr') then
+ nerr = nerr + 1
+ if (nerr .ge. maxerr) done = .true.
+ else
+ nerr = 0
+ end if
+!c
+!c test for too many total iterations
+!c
+ if (niter .ge. maxiter) then
+ statusbf = 'IterLimit'
+ done = .true.
+ end if
+!c
+!c test the normal termination criteria
+!c
+ if (f .le. fctmin) then
+ statusbf = 'SmallFct '
+ done = .true.
+ end if
+ if (g_rms .le. grdmin) then
+ statusbf = 'SmallGrad'
+ done = .true.
+ end if
+!c
+!c print intermediate results for the current iteration
+!c
+! write(iout,*) "niter1",niter
+ itime_mat=niter
+ if (jprint .gt. 0) then
+ if (done .or. mod(niter,jprint).eq.0) then
+ if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.&
+ g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.&
+ f_move.gt.-1.0d5) then
+ write (jout,80) niter,f,g_rms,f_move,x_move,&
+ angle,ncalls,statusbf
+ 80 format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
+ else
+ write (jout,90) niter,f,g_rms,f_move,x_move,&
+ angle,ncalls,statusbf
+ 90 format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
+ end if
+ end if
+ flush (jout)
+ end if
+!c
+!c write intermediate results for the current iteration
+!c
+ if (iwrite .gt. 0) then
+ if (done .or. mod(niter,iwrite).eq.0) then
+ call optsave (niter,f,x0)
+ end if
+ end if
+ end do
+!c
+!c perform deallocation of some local arrays
+!c
+ deallocate (x_old)
+ deallocate (g)
+ deallocate (g_old)
+ deallocate (p)
+ deallocate (q)
+ deallocate (r)
+ deallocate (h0)
+ if (msav .ne. 0) then
+ deallocate (rho)
+ deallocate (alpha)
+ deallocate (s)
+ deallocate (y)
+ end if
+!c
+!c set final value of the objective function
+!c
+ minimum = f
+ if (jprint .gt. 0) then
+ if (statusbf.eq.'SmallGrad' .or. statusbf.eq.'SmallFct ') then
+ write (jout,100) statusbf
+ 100 format (/,' LBFGS -- Normal Termination due to ',a9)
+ else
+ write (jout,110) statusbf
+ 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)
+ end if
+ flush (jout)
+ end if
+ return
+ end subroutine
+!ifdef LBFGS
+! double precision function funcgrad(x,g)
+! implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.MD'
+! include 'COMMON.QRESTR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! double precision energia(0:n_ene)
+! double precision x(nvar),g(nvar)
+! integer i
+! if (jjj.gt.0) then
+! write (iout,*) "in func x"
+! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+! endif
+! call var_to_geom(nvar,x)
+! call zerograd
+! call chainbuild_extconf
+! call etotal(energia(0))
+! call sum_gradient
+! funcgrad=energia(0)
+! call cart2intgrad(nvar,g)
+!
+! Add the components corresponding to local energy terms.
+!
+! Add the usampl contributions
+! if (usampl) then
+! do i=1,nres-3
+! gloc(i,icg)=gloc(i,icg)+dugamma(i)
+! enddo
+! do i=1,nres-2
+! gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+! enddo
+! endif
+! do i=1,nvar
+!d write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
+! g(i)=g(i)+gloc(i,icg)
+! enddo
+! return
+! end
+!endif