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. * * 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' integer icall common /srutu/ icall #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. 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 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) 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_extconf 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 #endif return end #ifdef MPI c---------------------------------------------------------------------------- subroutine ergastulum implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" double precision time00 integer ierr,ierror #endif include 'COMMON.SETUP' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.IOUNITS' 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) c double precision Gcopy(maxres2,maxres2) c 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) c write (*,*) 'Processor',fg_rank,' CG group',kolor, c & ' receives order from Master' #ifdef MPI 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 #endif 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 #ifndef FIVEDIAG 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) #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 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 ************************************************************************ #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' 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,*) "in func x" 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_extconf 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,*) "upon exit from func" c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n) c write (iout,*) 'f=',f 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_extconf 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------------------------------------------------------- #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 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 .and. itype(i).ne.ntyp1) then igall=igall+1 if (mask_side(i).eq.1) then ig=ig+1 xx(ig)=x(igall) 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 return end subroutine xx2x(x,xx) 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) 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 .and. itype(i).ne.ntyp1) then igall=igall+1 if (mask_side(i).eq.1) then ig=ig+1 x(igall)=xx(ig) 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 end 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 include 'COMMON.CONTROL' include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.VAR' 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 v(1:lv) common /przechowalnia/ v external func_dc,grad_dc,fdum #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 * 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 #endif 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 nvarx=k call chainbuild_cart 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 c aincr=1.0d-7 c call exec_checkgrad c----- call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum) c----- c write (iout,*) "checkgrad after SUMSL" c call exec_checkgrad c----- #endif 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 #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' #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 #endif