X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fminimize_p.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fminimize_p.F;h=06c7a7366ff90350b84f3948668b036b094ff1d7;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/minimize_p.F b/source/unres/src_MD-M-newcorr/minimize_p.F
new file mode 100644
index 0000000..06c7a73
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/minimize_p.F
@@ -0,0 +1,641 @@
+      subroutine minimize(etot,x,iretcode,nfun)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+*********************************************************************
+* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
+* the calling subprogram.                                           *     
+* when d(i)=1.0, then v(35) is the length of the initial step,      *     
+* calculated in the usual pythagorean way.                          *     
+* absolute convergence occurs when the function is within v(31) of  *     
+* zero. unless you know the minimum value in advance, abs convg     *     
+* is probably not useful.                                           *     
+* relative convergence is when the model predicts that the function *   
+* will decrease by less than v(32)*abs(fun).                        *   
+*********************************************************************
+      include 'COMMON.IOUNITS'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.MINIM'
+      common /srutu/ icall
+      dimension iv(liv)                                               
+      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+      double precision energia(0:n_ene)
+      external func,gradient,fdum
+      external func_restr,grad_restr
+      logical not_done,change,reduce 
+      common /przechowalnia/ v
+
+      icall = 1
+
+      NOT_DONE=.TRUE.
+
+c     DO WHILE (NOT_DONE)
+
+      call deflt(2,iv,liv,lv,v)                                         
+* 12 means fresh start, dont call deflt                                 
+      iv(1)=12                                                          
+* max num of fun calls                                                  
+      if (maxfun.eq.0) maxfun=500
+      iv(17)=maxfun
+* max num of iterations                                                 
+      if (maxmin.eq.0) maxmin=1000
+      iv(18)=maxmin
+* controls output                                                       
+      iv(19)=2                                                          
+* selects output unit                                                   
+      iv(21)=0
+      if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
+* 1 means to print out result                                           
+      iv(22)=print_min_res
+* 1 means to print out summary stats                                    
+      iv(23)=print_min_stat
+* 1 means to print initial x and d                                      
+      iv(24)=print_min_ini
+* min val for v(radfac) default is 0.1                                  
+      v(24)=0.1D0                                                       
+* max val for v(radfac) default is 4.0                                  
+      v(25)=2.0D0                                                       
+c     v(25)=4.0D0                                                       
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
+* the sumsl default is 0.1                                              
+      v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)                         
+* the sumsl default is 100*machep                                       
+      v(34)=v(34)/100.0D0                                               
+* absolute convergence                                                  
+      if (tolf.eq.0.0D0) tolf=1.0D-4
+      v(31)=tolf
+* relative convergence                                                  
+      if (rtolf.eq.0.0D0) rtolf=1.0D-4
+      v(32)=rtolf
+* controls initial step size                                            
+       v(35)=1.0D-1                                                    
+* large vals of d correspond to small components of step                
+      do i=1,nphi
+        d(i)=1.0D-1
+      enddo
+      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))
+c     etot = energia(0)
+      IF (mask_r) THEN
+       call x2xx(x,xx,nvar_restr)
+       call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
+     &                    iv,liv,lv,v,idum,rdum,fdum)      
+       call xx2x(x,xx)
+      ELSE
+       call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
+      ENDIF
+      etot=v(10)                                                      
+      iretcode=iv(1)
+cd    print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
+cd    write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
+c     call intout
+c     change=reduce(x)
+      call var_to_geom(nvar,x)
+c     if (change) then
+c       write (iout,'(a)') 'Reduction worked, minimizing again...'
+c     else
+c       not_done=.false.
+c     endif
+      call chainbuild
+c     call etotal(energia(0))
+c     etot=energia(0)
+c     call enerprint(energia(0))
+      nfun=iv(6)
+
+c     write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
+
+c     ENDDO ! NOT_DONE
+
+      return  
+      end  
+#ifdef MPI
+c----------------------------------------------------------------------------
+      subroutine ergastulum
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.FFIELD'
+      include 'COMMON.INTERACT'
+      include 'COMMON.MD'
+      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/
+C Workers wait for variables and NF, and NFL from the boss 
+      iorder=0
+      do while (iorder.ge.0)
+c      write (*,*) 'Processor',fg_rank,' CG group',kolor,
+c     & ' receives order from Master'
+        time00=MPI_Wtime()
+        call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
+        time_Bcast=time_Bcast+MPI_Wtime()-time00
+        if (icall.gt.4 .and. iorder.ge.0) 
+     &   time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
+       icall=icall+1
+c      write (*,*) 
+c     & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
+        if (iorder.eq.0) then
+          call zerograd
+          call etotal(edum)
+c          write (2,*) "After etotal"
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+c          call flush(2)
+        else if (iorder.eq.2) then
+          call zerograd
+          call etotal_short(edum)
+c          write (2,*) "After etotal_short"
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+c          call flush(2)
+        else if (iorder.eq.3) then
+          call zerograd
+          call etotal_long(edum)
+c          write (2,*) "After etotal_long"
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+c          call flush(2)
+        else if (iorder.eq.1) then
+          call sum_gradient
+c          write (2,*) "After sum_gradient"
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+c          call flush(2)
+        else if (iorder.eq.4) then
+          call ginv_mult(z,d_a_tmp)
+        else if (iorder.eq.5) then
+c Setup MD things for a slave
+          dimen=(nct-nnt+1)+nside
+          dimen1=(nct-nnt)+(nct-nnt+1)
+          dimen3=dimen*3
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+c          call flush(2)
+          call int_bounds(dimen,igmult_start,igmult_end)
+          igmult_start=igmult_start-1
+          call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
+     &     ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
+           my_ng_count=igmult_end-igmult_start
+          call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
+     &     MPI_INTEGER,FG_COMM,IERROR)
+c          write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
+c          write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
+          myginv_ng_count=maxres2*my_ng_count
+c          write (2,*) "igmult_start",igmult_start," igmult_end",
+c     &     igmult_end," my_ng_count",my_ng_count
+c          call flush(2)
+          call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
+     &     nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
+          call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
+     &     nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
+c          write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
+c          write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
+c          call flush(2)
+c          call MPI_Barrier(FG_COMM,IERROR)
+          time00=MPI_Wtime()
+          call MPI_Scatterv(ginv(1,1),nginv_counts(0),
+     &     nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
+     &     myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
+#ifdef TIMING
+          time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
+#endif
+          do i=1,dimen
+            do j=1,2*my_ng_count
+              ginv(j,i)=gcopy(i,j)
+            enddo
+          enddo
+c          write (2,*) "dimen",dimen," dimen3",dimen3
+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)
+        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
+        else if (iorder.eq.9) then
+          call fricmat_mult(z,d_a_tmp)
+        else if (iorder.eq.10) then
+          call setup_fricmat
+        endif
+      enddo
+      write (*,*) 'Processor',fg_rank,' CG group',kolor,
+     &  ' absolute rank',myrank,' leves ERGASTULUM.'
+      write(*,*)'Processor',fg_rank,' wait times for respective orders',
+     & (' order[',i,']',time_order(i),i=0,10)
+      return
+      end
+#endif
+************************************************************************
+      subroutine func(n,x,nf,f,uiparm,urparm,ufparm)  
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      common /chuju/ jjj
+      double precision energia(0:n_ene)
+      integer jjj
+      double precision ufparm
+      external ufparm                                                   
+      integer uiparm(1)                                                 
+      real*8 urparm(1)                                                    
+      dimension x(maxvar)
+c     if (jjj.gt.0) then
+c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c     endif
+      nfl=nf
+      icg=mod(nf,2)+1
+cd      print *,'func',nf,nfl,icg
+      call var_to_geom(n,x)
+      call zerograd
+      call chainbuild
+cd    write (iout,*) 'ETOTAL called from FUNC'
+      call etotal(energia(0))
+      call sum_gradient
+      f=energia(0)
+c     if (jjj.gt.0) then
+c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c       write (iout,*) 'f=',etot
+c       jjj=0
+c     endif
+      return                                                            
+      end                                                               
+************************************************************************
+      subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)  
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      common /chuju/ jjj
+      double precision energia(0:n_ene)
+      integer jjj
+      double precision ufparm
+      external ufparm                                                   
+      integer uiparm(1)                                                 
+      real*8 urparm(1)                                                    
+      dimension x(maxvar)
+c     if (jjj.gt.0) then
+c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c     endif
+      nfl=nf
+      icg=mod(nf,2)+1
+      call var_to_geom_restr(n,x)
+      call zerograd
+      call chainbuild
+cd    write (iout,*) 'ETOTAL called from FUNC'
+      call etotal(energia(0))
+      call sum_gradient
+      f=energia(0)
+c     if (jjj.gt.0) then
+c       write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c       write (iout,*) 'f=',etot
+c       jjj=0
+c     endif
+      return                                                            
+      end                                                               
+c-------------------------------------------------------
+      subroutine x2xx(x,xx,n)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      double precision xx(maxvar),x(maxvar)
+
+      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                                                               
+          xx(ig)=x(igall)                       
+        endif                                                                   
+      enddo                                                                     
+                                                                                
+      do i=3,nres                                                               
+        igall=igall+1                                                           
+        if (mask_theta(i).eq.1) then                                            
+          ig=ig+1                                                               
+          xx(ig)=x(igall)
+        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                                                             
+            xx(ig)=x(igall)
+          endif                                                                 
+        endif                                                                   
+      enddo                                                                     
+      enddo                              
+ 
+      n=ig
+
+      return
+      end
+
+      subroutine xx2x(x,xx)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      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                                                               
+          x(igall)=xx(ig)
+        endif                                                                   
+      enddo                                                                     
+                                                                                
+      do i=3,nres                                                               
+        igall=igall+1                                                           
+        if (mask_theta(i).eq.1) then                                            
+          ig=ig+1                                                               
+          x(igall)=xx(ig)
+        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                                                             
+            x(igall)=xx(ig)
+          endif                                                                 
+        endif                                                                   
+      enddo                                                             
+      enddo                              
+
+      return
+      end
+  
+c---------------------------------------------------------- 
+      subroutine minim_dc(etot,iretcode,nfun)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.MINIM'
+      include 'COMMON.CHAIN'
+      dimension iv(liv)                                               
+      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+      common /przechowalnia/ v
+
+      double precision energia(0:n_ene)
+      external func_dc,grad_dc,fdum
+      logical not_done,change,reduce 
+      double precision g(maxvar),f1
+
+      call deflt(2,iv,liv,lv,v)                                         
+* 12 means fresh start, dont call deflt                                 
+      iv(1)=12                                                          
+* max num of fun calls                                                  
+      if (maxfun.eq.0) maxfun=500
+      iv(17)=maxfun
+* max num of iterations                                                 
+      if (maxmin.eq.0) maxmin=1000
+      iv(18)=maxmin
+* controls output                                                       
+      iv(19)=2                                                          
+* selects output unit                                                   
+      iv(21)=0
+      if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
+* 1 means to print out result                                           
+      iv(22)=print_min_res
+* 1 means to print out summary stats                                    
+      iv(23)=print_min_stat
+* 1 means to print initial x and d                                      
+      iv(24)=print_min_ini
+* min val for v(radfac) default is 0.1                                  
+      v(24)=0.1D0                                                       
+* max val for v(radfac) default is 4.0                                  
+      v(25)=2.0D0                                                       
+c     v(25)=4.0D0                                                       
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
+* the sumsl default is 0.1                                              
+      v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)                         
+* the sumsl default is 100*machep                                       
+      v(34)=v(34)/100.0D0                                               
+* absolute convergence                                                  
+      if (tolf.eq.0.0D0) tolf=1.0D-4
+      v(31)=tolf
+* relative convergence                                                  
+      if (rtolf.eq.0.0D0) rtolf=1.0D-4
+      v(32)=rtolf
+* controls initial step size                                            
+       v(35)=1.0D-1                                                    
+* large vals of d correspond to small components of step                
+      do i=1,6*nres
+        d(i)=1.0D-1
+      enddo
+
+      k=0
+      do i=1,nres-1
+        do j=1,3
+          k=k+1
+          x(k)=dc(j,i)
+        enddo
+      enddo
+      do i=2,nres-1
+        if (ialph(i,1).gt.0) then
+        do j=1,3
+          k=k+1
+          x(k)=dc(j,i+nres)
+        enddo
+        endif
+      enddo
+
+      call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)      
+
+      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
+
+cd      call zerograd
+cd      nf=0
+cd      call func_dc(k,x,nf,f,idum,rdum,fdum)
+cd      call grad_dc(k,x,nf,g,idum,rdum,fdum)
+cd
+cd      do i=1,k
+cd       x(i)=x(i)+1.0D-5
+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
+
+      etot=v(10)                                                      
+      iretcode=iv(1)
+      nfun=iv(6)
+      return  
+      end  
+
+      subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)  
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      double precision energia(0:n_ene)
+      double precision ufparm
+      external ufparm                                                   
+      integer uiparm(1)                                                 
+      real*8 urparm(1)                                                    
+      dimension x(maxvar)
+      nfl=nf
+cbad      icg=mod(nf,2)+1
+      icg=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))
+      f=energia(0)
+
+cd      print *,'func_dc ',nf,nfl,f
+
+      return                                                            
+      end                                                               
+
+      subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
+      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'
+      external ufparm
+      integer uiparm(1),k
+      double precision urparm(1)
+      dimension x(maxvar),g(maxvar)
+c
+c
+c
+cbad      icg=mod(nf,2)+1
+      icg=1
+cd      print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
+      if (nf-nfl+1) 20,30,40
+   20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
+cd      print *,20
+      if (nf.eq.0) return
+      goto 40
+   30 continue
+cd      print *,30
+      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
+
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+   40 call cartgrad
+cd      print *,40
+      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