X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fminimize_p.F;h=a56e4f88de78b514352c594bb189b373642dc8c6;hb=a30bd29e64da2aa47b84963fdd0bf4192ead2738;hp=f9faf7c2bf08e37cef6b326257f18868e6801ba0;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;p=unres.git diff --git a/source/unres/src-HCD-5D/minimize_p.F b/source/unres/src-HCD-5D/minimize_p.F index f9faf7c..a56e4f8 100644 --- a/source/unres/src-HCD-5D/minimize_p.F +++ b/source/unres/src-HCD-5D/minimize_p.F @@ -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. * @@ -17,15 +27,61 @@ 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