make cp src-HCD-5D
[unres.git] / source / unres / src-HCD-5D / minimize_p.F
index f9faf7c..a56e4f8 100644 (file)
@@ -1,7 +1,17 @@
       subroutine minimize(etot,x,iretcode,nfun)
-      implicit real*8 (a-h,o-z)
+#ifdef LBFGS
+      use minima
+      use inform
+      use output
+      use iounit
+      use scales
+#endif
+      implicit none
       include 'DIMENSIONS'
+#ifndef LBFGS
+      integer liv,lv
       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+#endif
 *********************************************************************
 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
 * the calling subprogram.                                           *     
       include 'COMMON.VAR'
       include 'COMMON.GEO'
       include 'COMMON.MINIM'
+      integer icall
       common /srutu/ icall
-      dimension iv(liv)                                               
-      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
-      double precision energia(0:n_ene)
+#ifdef LBFGS
+      double precision grdmin
+      external funcgrad
+      external optsave
+#else
+      integer iv(liv)                                               
+      double precision v(1:lv)
+      common /przechowalnia/ v
+      integer idum
+      double precision rdum
+      double precision fdum
       external func,gradient,fdum
       external func_restr,grad_restr
       logical not_done,change,reduce 
+#endif
+      double precision x(maxvar),d(maxvar),xx(maxvar)
+      double precision energia(0:n_ene)
+      integer i,nvar_restr,nfun,iretcode
+      double precision etot
 c      common /przechowalnia/ v
 
+#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
+c      write (iout,*) "Calling lbfgs"
+      write (iout,*) 'Calling LBFGS, minimization in angles'
+      call var_to_geom(nvar,x)
+      call chainbuild_extconf
+      call etotal(energia(0))
+      call enerprint(energia(0))
+      call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
+      deallocate(scale)
+      write (iout,*) "Minimized energy",etot
+#else
       icall = 1
 
       NOT_DONE=.TRUE.
@@ -78,10 +134,12 @@ c     v(25)=4.0D0
       do i=nphi+1,nvar
         d(i)=1.0D-1
       enddo
-cd    print *,'Calling SUMSL'
-c     call var_to_geom(nvar,x)
-c     call chainbuild
-c     call etotal(energia(0))
+      write (iout,*) 'Calling SUMSL'
+      call var_to_geom(nvar,x)
+      call chainbuild_extconf
+      call intout
+      call etotal(energia(0))
+      call enerprint(energia(0))
 c     etot = energia(0)
       IF (mask_r) THEN
        call x2xx(x,xx,nvar_restr)
@@ -103,7 +161,7 @@ c       write (iout,'(a)') 'Reduction worked, minimizing again...'
 c     else
 c       not_done=.false.
 c     endif
-      call chainbuild
+      call chainbuild_extconf
 c     call etotal(energia(0))
 c     etot=energia(0)
 c     call enerprint(energia(0))
@@ -112,16 +170,18 @@ c     call enerprint(energia(0))
 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
 
 c     ENDDO ! NOT_DONE
-
+#endif
       return  
       end  
 #ifdef MPI
 c----------------------------------------------------------------------------
       subroutine ergastulum
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include "mpif.h"
+      double precision time00
+      integer ierr,ierror
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.DERIV'
@@ -130,12 +190,18 @@ c----------------------------------------------------------------------------
       include 'COMMON.FFIELD'
       include 'COMMON.INTERACT'
       include 'COMMON.MD'
+#ifdef FIVEDIAG
+       include 'COMMON.LAGRANGE.5diag'
+#else
+       include 'COMMON.LAGRANGE'
+#endif
       include 'COMMON.TIME1'
       double precision z(maxres6),d_a_tmp(maxres6)
       double precision edum(0:n_ene),time_order(0:10)
       double precision Gcopy(maxres2,maxres2)
       common /przechowalnia/ Gcopy
       integer icall /0/
+      integer i,j,iorder
 C Workers wait for variables and NF, and NFL from the boss 
       iorder=0
       do while (iorder.ge.0)
@@ -173,7 +239,8 @@ c          call flush(2)
           call sum_gradient
 c          write (2,*) "After sum_gradient"
 c          write (2,*) "dimen",dimen," dimen3",dimen3
-c          call flush(2)
+c          call flush(2
+#ifndef FIVEDIAG
         else if (iorder.eq.4) then
           call ginv_mult(z,d_a_tmp)
         else if (iorder.eq.5) then
@@ -221,14 +288,17 @@ c          write (2,*) "End MD setup"
 c          call flush(2)
 c           write (iout,*) "My chunk of ginv_block"
 c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
+#endif
         else if (iorder.eq.6) then
           call int_from_cart1(.false.)
         else if (iorder.eq.7) then
           call chainbuild_cart
         else if (iorder.eq.8) then
           call intcartderiv
+#ifndef FIVEDIAG
         else if (iorder.eq.9) then
           call fricmat_mult(z,d_a_tmp)
+#endif
         else if (iorder.eq.10) then
           call setup_fricmat
         endif
@@ -241,6 +311,53 @@ c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
       end
 #endif
 ************************************************************************
+#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
+c     if (jjj.gt.0) then
+c      write (iout,*) "in func x"
+c      write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c     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)
+C
+C Add the components corresponding to local energy terms.
+C
+c 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
+cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
+        g(i)=g(i)+gloc(i,icg)
+      enddo
+      return                                                            
+      end                                                               
+#else
       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)  
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -312,46 +429,51 @@ c     endif
       return                                                            
       end                                                               
 c-------------------------------------------------------
+#endif
       subroutine x2xx(x,xx,n)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
       double precision xx(maxvar),x(maxvar)
 
+c      write (iout,*) "nvar",nvar
       do i=1,nvar
         varall(i)=x(i)
       enddo
 
-      ig=0                                                                      
-      igall=0                                                                   
-      do i=4,nres                                                               
-        igall=igall+1                                                           
-        if (mask_phi(i).eq.1) then                                              
-          ig=ig+1                                                               
+      ig=0                                             
+      igall=0                                          
+      do i=4,nres                                      
+        igall=igall+1                                  
+        if (mask_phi(i).eq.1) then                     
+          ig=ig+1                                      
           xx(ig)=x(igall)                       
-        endif                                                                   
-      enddo                                                                     
-                                                                                
-      do i=3,nres                                                               
-        igall=igall+1                                                           
-        if (mask_theta(i).eq.1) then                                            
-          ig=ig+1                                                               
+        endif                                          
+      enddo                                            
+                                                       
+      do i=3,nres                                      
+        igall=igall+1                                  
+        if (mask_theta(i).eq.1) then                   
+          ig=ig+1                                      
           xx(ig)=x(igall)
-        endif                                                                   
+        endif                                          
       enddo                                          
 
-      do ij=1,2                                                                 
-      do i=2,nres-1                                                             
-        if (itype(i).ne.10) then                                                
-          igall=igall+1                                                         
-          if (mask_side(i).eq.1) then                                           
-            ig=ig+1                                                             
+      do ij=1,2                                        
+      do i=2,nres-1                                    
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+          igall=igall+1                                 
+          if (mask_side(i).eq.1) then                   
+            ig=ig+1                                     
             xx(ig)=x(igall)
-          endif                                                                 
-        endif                                                                   
-      enddo                                                                     
+c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
+c            write (iout,*) "x",x(igall)," xx",xx(ig)
+          endif                                         
+        endif                                           
+      enddo                                             
       enddo                              
  
       n=ig
@@ -365,40 +487,43 @@ c-------------------------------------------------------
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
       double precision xx(maxvar),x(maxvar)
 
       do i=1,nvar
         x(i)=varall(i)
       enddo
 
-      ig=0                                                                      
-      igall=0                                                                   
-      do i=4,nres                                                               
-        igall=igall+1                                                           
-        if (mask_phi(i).eq.1) then                                              
-          ig=ig+1                                                               
+      ig=0                                                     
+      igall=0                                                  
+      do i=4,nres                                              
+        igall=igall+1                                          
+        if (mask_phi(i).eq.1) then                             
+          ig=ig+1                                              
           x(igall)=xx(ig)
-        endif                                                                   
-      enddo                                                                     
-                                                                                
-      do i=3,nres                                                               
-        igall=igall+1                                                           
-        if (mask_theta(i).eq.1) then                                            
-          ig=ig+1                                                               
+        endif                                                  
+      enddo                                                    
+                                                               
+      do i=3,nres                                              
+        igall=igall+1                                          
+        if (mask_theta(i).eq.1) then                           
+          ig=ig+1                                              
           x(igall)=xx(ig)
-        endif                                                                   
+        endif                                                  
       enddo                                          
 
-      do ij=1,2                                                                 
-      do i=2,nres-1                                                             
-        if (itype(i).ne.10) then                                                
-          igall=igall+1                                                         
-          if (mask_side(i).eq.1) then                                           
-            ig=ig+1                                                             
+      do ij=1,2                                                
+      do i=2,nres-1                                            
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+          igall=igall+1                                        
+          if (mask_side(i).eq.1) then                          
+            ig=ig+1                                            
             x(igall)=xx(ig)
-          endif                                                                 
-        endif                                                                   
-      enddo                                                             
+c            write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
+c            write (iout,*) "x",x(igall)," xx",xx(ig)
+          endif                                                
+        endif                                                  
+      enddo                                             
       enddo                              
 
       return
@@ -406,9 +531,18 @@ c-------------------------------------------------------
   
 c---------------------------------------------------------- 
       subroutine minim_dc(etot,iretcode,nfun)
+#ifdef LBFGS
+      use minima
+      use inform
+      use output
+      use iounit
+      use scales
+#endif
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+#ifndef LBFGS
       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+#endif
 #ifdef MPI
       include 'mpif.h'
 #endif
@@ -419,15 +553,28 @@ c----------------------------------------------------------
       include 'COMMON.GEO'
       include 'COMMON.MINIM'
       include 'COMMON.CHAIN'
+      double precision minval,x(maxvar),d(maxvar),xx(maxvar)
+#ifdef LBFGS
+      double precision grdmin
+      double precision funcgrad_dc
+      external funcgrad_dc,optsave
+#else
       dimension iv(liv)                                               
-      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+      double precision v(1:lv)
       common /przechowalnia/ v
-
-      double precision energia(0:n_ene)
       external func_dc,grad_dc,fdum
-      logical not_done,change,reduce 
+#endif
       double precision g(maxvar),f1
-
+      integer nvarx
+      double precision energia(0:n_ene)
+#ifdef LBFGS
+      maxiter=maxmin
+      coordtype='CARTESIAN'
+      grdmin=tolf
+      jout=iout
+      jprint=print_min_stat
+      iwrite=0
+#else
       call deflt(2,iv,liv,lv,v)                                         
 * 12 means fresh start, dont call deflt                                 
       iv(1)=12                                                          
@@ -471,7 +618,7 @@ c     v(25)=4.0D0
       do i=1,6*nres
         d(i)=1.0D-1
       enddo
-
+#endif
       k=0
       do i=1,nres-1
         do j=1,3
@@ -487,6 +634,37 @@ c     v(25)=4.0D0
         enddo
         endif
       enddo
+      nvarx=k
+      write (iout,*) "Variables set up nvarx",nvarx
+      write (iout,*) "Before energy minimization"
+      call etotal(energia(0))
+      call enerprint(energia(0))
+#ifdef LBFGS
+c
+c     From tinker
+c
+c     perform dynamic allocation of some global arrays
+c
+      if (.not. allocated(scale))  allocate (scale(nvarx))
+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, nvarx
+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
+c      write (iout,*) "minim_dc Calling lbfgs"
+      call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
+      deallocate(scale)
+c      write (iout,*) "minim_dc After lbfgs"
+#else
 c-----
 c      write (iout,*) "checkgrad before SUMSL"
 c      icheckgrad=1
@@ -499,7 +677,7 @@ c-----
 c      write (iout,*) "checkgrad after SUMSL"
 c      call exec_checkgrad
 c-----
-
+#endif
       k=0
       do i=1,nres-1
         do j=1,3
@@ -528,13 +706,77 @@ cd       call func_dc(k,x,nf,f1,idum,rdum,fdum)
 cd       x(i)=x(i)-1.0D-5
 cd       print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
 cd      enddo
-
+#ifndef LBFGS
       etot=v(10)                                                      
       iretcode=iv(1)
       nfun=iv(6)
+#endif
       return  
       end  
+#ifdef LBFGS
+      double precision function funcgrad_dc(x,g)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.MD'
+      include 'COMMON.IOUNITS'
+      integer k
+      dimension x(maxvar),g(maxvar)
+      double precision energia(0:n_ene)
+      common /gacia/ nf
+c
+      nf=nf+1
+      k=0
+      do i=1,nres-1
+        do j=1,3
+          k=k+1
+          dc(j,i)=x(k)
+        enddo
+      enddo
+      do i=2,nres-1
+        if (ialph(i,1).gt.0) then
+        do j=1,3
+          k=k+1
+          dc(j,i+nres)=x(k)
+        enddo
+        endif
+      enddo
+      call chainbuild_cart
+      call zerograd
+      call etotal(energia(0))
+c      write (iout,*) "energia",energia(0)
+      funcgrad_dc=energia(0)
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+      call cartgrad
+      k=0
+      do i=1,nres-1
+        do j=1,3
+          k=k+1
+          g(k)=gcart(j,i)
+        enddo
+      enddo
+      do i=2,nres-1
+        if (ialph(i,1).gt.0) then
+        do j=1,3
+          k=k+1
+          g(k)=gxcart(j,i)
+        enddo
+        endif
+      enddo       
 
+      return
+      end
+#else
       subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)  
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -551,7 +793,7 @@ cd      enddo
       double precision ufparm
       external ufparm                                                   
       integer uiparm(1)                                                 
-      real*8 urparm(1)                                                    
+      real*8 urparm(1)
       dimension x(maxvar)
       nfl=nf
 cbad      icg=mod(nf,2)+1
@@ -654,3 +896,4 @@ cd      print *,40
 
       return
       end
+#endif