make cp src-HCD-5D

[unres.git] / source / unres / src-HCD-5D / minim_jlee.F

diff --git a/source/unres/src-HCD-5D/minim_jlee.F b/source/unres/src-HCD-5D/minim_jlee.F
index 56d5010..7162afb 100644 (file)
--- a/source/unres/src-HCD-5D/minim_jlee.F
+++ b/source/unres/src-HCD-5D/minim_jlee.F
@@ -1,13 +1,33 @@
       subroutine minim_jlee
+#ifdef LBFGS
+      use minima
+      use inform
+      use output
+      use iounit
+      use scales
+#endif
 c  controls minimization and sorting routines
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
-      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+#ifndef LBFGS
+      integer liv,lv
+      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+#endif
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       include 'COMMON.MINIM'
       include 'COMMON.CONTROL'
+#ifdef LBFGS
+      common /gacia/ nfun
+      double precision grdmin
+      external funcgrad
+      external optsave
+#else
       external func,gradient,fdum
+      dimension iv(liv)                                               
+      double precision v(1:lv+1)
+      common /przechowalnia/ v
+#endif
       real ran1,ran2,ran3
 #ifdef MPI
       include 'mpif.h'
@@ -22,19 +42,40 @@ c  controls minimization and sorting routines
       dimension var(maxvar),erg(mxch*(mxch+1)/2+1)
       dimension var2(maxvar)
       integer iffr(maxres),ihpbt(maxdim),jhpbt(maxdim)
-      double precision d(maxvar),v(1:lv+1),garbage(maxvar)                     
+      double precision d(maxvar),garbage(maxvar),g(maxvar)
       double precision energia(0:n_ene),time0s,time1s
       dimension indx(9),info(12)
-      dimension iv(liv)                                               
       dimension idum(1),rdum(1)
       dimension icont(2,maxcont)
       logical check_var,fail
       integer iloop(2)
-      common /przechowalnia/ v
       data rad /1.745329252d-2/
 c  receive # of start
 !      print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
 !     &   ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
+#ifdef LBFGS
+      maxiter=maxmin
+      coordtype='RIGIDBODY'
+      grdmin=tolf
+      jout=iout
+      jprint=print_min_stat
+      iwrite=0
+      if (.not. allocated(scale))  allocate (scale(nvar))
+c
+c     set scaling parameter for function and derivative values;
+c     use square root of median eigenvalue of typical Hessian
+c
+      set_scale = .true.
+c      nvar = 0
+      do i = 1, nvar
+c         if (use(i)) then
+c            do j = 1, 3
+c               nvar = nvar + 1
+               scale(i) = 12.0d0
+c            end do
+c         end if
+      end do
+#endif
       nhpb0=nhpb
    10 continue
       time0s=MPI_WTIME()
@@ -161,8 +202,13 @@ crc overlap test
              nfun=nfun+1
              write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
             else
+#ifdef LBFGS
+             etot=1.0d20
+             nfun=-1
+#else
              v(10)=1.0d20
              iv(1)=-1
+#endif
              goto 201
             endif
           endif
@@ -176,8 +222,12 @@ cd          write(iout,*) 'sc_move',nft_sc,etot
       endif 
 
       if (check_var(var,info)) then 
+#ifdef LBFGS
+           etot=1.0d21
+#else
            v(10)=1.0d21
            iv(1)=6
+#endif
            goto 201
       endif
 
@@ -189,10 +239,22 @@ crc
 !      write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
 !      write (*,'(8f10.4)') (var(i),i=1,nvar)
 
-       do i=1,nvar
-         garbage(i)=var(i)
-       enddo
+      do i=1,nvar
+        garbage(i)=var(i)
+      enddo
+#ifdef LBFGS
+      eee=funcgrad(var,g)
+      nfun=nfun+1
+      if(eee.ge.1.0d20) then
+c       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
+c       print *,' energy before SUMSL =',eee
+c       print *,' aborting local minimization'
+       go to 201
+      endif
 
+      call lbfgs (nvar,var,etot,grdmin,funcgrad,optsave)
+      deallocate(scale)
+#else
       call deflt(2,iv,liv,lv,v)                                         
 * 12 means fresh start, dont call deflt                                 
       iv(1)=12                                                          
@@ -262,8 +324,12 @@ c      print *, 'MINIM_JLEE: ',me,' before SUMSL '
 c       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
 c       print *,' energy before SUMSL =',eee
 c       print *,' aborting local minimization'
+#ifdef LBFGS
+       etot=eee
+#else
        iv(1)=-1
        v(10)=eee
+#endif
        go to 201
       endif
 
@@ -274,6 +340,7 @@ c      print *, 'MINIM_JLEE: ',me,' after SUMSL '
 
 c  find which conformation was returned from sumsl
         nfun=nfun+iv(7)
+#endif
 !      print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
 !     & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
 c        if (iv(1).ne.4 .or. nf.le.1) then
@@ -311,7 +378,11 @@ c       print *, 'MINIM_JLEE: ',me,' minimized: ',n
   201  continue
         indx(1)=n
 c return code: 6-gradient 9-number of ftn evaluation, etc
+#ifdef LBFGS
+        indx(2)=nfun
+#else
         indx(2)=iv(1)
+#endif
 c total # of ftn evaluations (for iwf=0, it includes all minimizations).
         indx(3)=nfun
         indx(4)=info(2)
@@ -325,12 +396,21 @@ c total # of ftn evaluations (for iwf=0, it includes all minimizations).
 c  send back energies
 c al & cc
 c calculate contact order
+#ifdef LBFGS
+#ifdef CO_BIAS
+        call contact(.false.,ncont,icont,co)
+        erg(1)=etot-1.0d2*co
+#else
+        erg(1)=etot
+#endif
+#else
 #ifdef CO_BIAS
         call contact(.false.,ncont,icont,co)
         erg(1)=v(10)-1.0d2*co
 #else
         erg(1)=v(10)
 #endif
+#endif
         j=1
         call mpi_send(erg,j,mpi_double_precision,king,idreal,
      *                 CG_COMM,ierr)