adam's update
[unres.git] / source / unres / src-HCD-5D / minimize_p.F
index f163846..6b9d204 100644 (file)
@@ -1,8 +1,17 @@
       subroutine minimize(etot,x,iretcode,nfun)
+#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.MINIM'
       integer icall
       common /srutu/ icall
-      integer 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
+      dimension 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.
@@ -85,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)
@@ -110,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))
@@ -119,7 +170,7 @@ c     call enerprint(energia(0))
 c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
 
 c     ENDDO ! NOT_DONE
-
+#endif
       return  
       end  
 #ifdef MPI
@@ -260,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'
@@ -331,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
@@ -384,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
@@ -451,7 +557,7 @@ c----------------------------------------------------------
 #ifdef LBFGS
       double precision grdmin
       double precision funcgrad_dc
-      external funcgrad_dc
+      external funcgrad_dc,optsave
 #else
       dimension iv(liv)                                               
       double precision v(1:lv)
@@ -550,9 +656,10 @@ c               nvar = nvar + 1
 c            end do
 c         end if
       end do
-      write (iout,*) "Calling lbfgs"
+c      write (iout,*) "minim_dc Calling lbfgs"
       call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
-      write (iout,*) "After lbfgs"
+      deallocate(scale)
+c      write (iout,*) "minim_dc After lbfgs"
 #else
 c-----
 c      write (iout,*) "checkgrad before SUMSL"
@@ -620,7 +727,9 @@ cd      enddo
       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
@@ -645,7 +754,6 @@ C
 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 C
       call cartgrad
-cd      print *,40
       k=0
       do i=1,nres-1
         do j=1,3