From: Adam Sieradzan Date: Tue, 13 Jun 2023 07:54:00 +0000 (+0200) Subject: cleaning water X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;ds=sidebyside;h=9591847f46e1b0163f508e2fda6a240fdc3e4890;p=unres4.git cleaning water --- diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt index 3fa4bc8..d9be5ea 100644 --- a/source/unres/CMakeLists.txt +++ b/source/unres/CMakeLists.txt @@ -22,11 +22,21 @@ set(UNRES_MD_SRC0 data/MPI_data.F90 data/REMD_data.F90 data/comm_local.F90 + inform.F90 + keys.F90 + math2.F90 + iounit.F90 + linmin.F90 + minima.F90 + output.F90 + scales.F90 + search.F90 + optsave_dum.F90 prng_32.F90 math.F90 random.F90 geometry.F90 - io_base.F90 + io_base.F90 energy.F90 check_bond.F90 control.F90 @@ -72,7 +82,7 @@ if (Fortran_COMPILER_NAME STREQUAL "ifort") set (CMAKE_Fortran_FLAGS_RELEASE " ") set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -traceback") # set(FFLAGS0 "-fpp -c -CB -g -ip " ) - set(FFLAGS0 "-O3 -g -ip -fpp -mcmodel=large" ) + set(FFLAGS0 "-CB -g -ip -fpp -mcmodel=large" ) # set(FFLAGS0 "-O0 -CB -CA -g" ) set(FFLAGS1 "-fpp -c -O " ) set(FFLAGS2 "-fpp -c -g -CA -CB ") @@ -125,7 +135,7 @@ endif(UNRES_MD_FF STREQUAL "GAB") if(UNRES_NEWGRAD) # set(CPPFLAGS "${CPPFLAGS} -DFIVEDIAG ") -set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG -DSC_END") +set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG -DSC_END -DLBFGS") endif() # set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG ") diff --git a/source/unres/CSA.F90 b/source/unres/CSA.F90 index 9e7dbfa..826f37e 100644 --- a/source/unres/CSA.F90 +++ b/source/unres/CSA.F90 @@ -70,6 +70,7 @@ !----------------------------------------------------------------------------- ! Maximum number of moves (n1-n8) integer,parameter :: mxmv=18 + !----------------------------------------------------------------------------- ! ! @@ -1028,7 +1029,14 @@ ! minim_jlee.F !----------------------------------------------------------------------------- subroutine minim_jlee - +#ifdef LBFGS + use minima + use inform + use output + use iounit + use scales + use energy, only: funcgrad +#endif use minim_data use MPI_data use energy_data @@ -1039,6 +1047,9 @@ use energy, only:fdum use control, only:init_int_table use minimm, only:sumsl,deflt +#ifdef LBFGS + use minimm,only: lbfgs +#endif ! controls minimization and sorting routines ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -1049,6 +1060,11 @@ include 'mpif.h' integer,parameter :: liv=60 integer :: lv +#ifdef LBFGS + double precision grdmin + real(kind=8),dimension(6*nres)::g + external optsave +#endif ! external func,gradient!,fdum !use minim & energy ! real(kind=4) :: ran1,ran2,ran3 ! include 'COMMON.SETUP' @@ -1075,7 +1091,10 @@ logical :: fail !check_var, integer :: iloop(2) !el common /przechowalnia/ v - integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode + integer :: i,j,ierr,n,nft_sc,nf,ierror,ierrcode +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: rad,eee,etot !,fdum !el from subroutine parmread ! Define the constants of the disulfide bridge @@ -1088,7 +1107,33 @@ ! receive # of start ! print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun, ! & ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf +#ifndef LBFGS if (.not. allocated(v)) allocate(v(1:lv)) +#endif +#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 +#endif + nhpb0=nhpb 10 continue time0s=MPI_WTIME() @@ -1215,8 +1260,15 @@ nfun=nfun+1 write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1) else +! v(10)=1.0d20 +! iv(1)=-1 +#ifdef LBFGS + etot=1.0d20 + nfun=-1 +#else v(10)=1.0d20 iv(1)=-1 +#endif goto 201 endif endif @@ -1230,8 +1282,15 @@ endif if (check_var(var,info)) then +! v(10)=1.0d21 +! iv(1)=6 +#ifdef LBFGS + etot=1.0d21 +#else v(10)=1.0d21 iv(1)=6 +#endif + goto 201 endif @@ -1246,6 +1305,19 @@ do i=1,nvar garbage(i)=var(i) enddo +#ifdef LBFGS + eee=funcgrad(var,g) + nfun=nfun+1 + if(eee.ge.1.0d20) then +!c print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL' +!c print *,' energy before SUMSL =',eee +!c print *,' aborting local minimization' + go to 201 + endif + + call lbfgs (nvar,var,etot,grdmin,funcgrad,optsave) + deallocate(scale) +#else call deflt(2,iv,liv,lv,v) ! 12 means fresh start, dont call deflt @@ -1328,6 +1400,7 @@ ! find which conformation was returned from sumsl nfun=nfun+iv(7) +#endif ! print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf, ! & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32) ! if (iv(1).ne.4 .or. nf.le.1) then @@ -1365,7 +1438,12 @@ 201 continue indx(1)=n ! return code: 6-gradient 9-number of ftn evaluation, etc +! indx(2)=iv(1) +#ifdef LBFGS + indx(2)=nfun +#else indx(2)=iv(1) +#endif ! total # of ftn evaluations (for iwf=0, it includes all minimizations). indx(3)=nfun indx(4)=info(2) @@ -1379,12 +1457,21 @@ ! send back energies ! al & cc ! calculate contact order +#ifdef LBFGS #ifdef CO_BIAS - call contact(.false.,ncont,icont_,co) + call contact(.false.,nnt,nct,ncont,icont,co) + erg(1)=etot-1.0d2*co +#else + erg(1)=etot +#endif +#else +#ifdef CO_BIAS + call contact(.false.,nnt,nct,ncont,icont,co) erg(1)=v(10)-1.0d2*co #else erg(1)=v(10) #endif +#endif j=1 call mpi_send(erg,j,mpi_double_precision,king,idreal,& CG_COMM,ierr) diff --git a/source/unres/MCM_MD.F90 b/source/unres/MCM_MD.F90 index eb07d69..090a7d6 100644 --- a/source/unres/MCM_MD.F90 +++ b/source/unres/MCM_MD.F90 @@ -344,9 +344,12 @@ real(kind=8),dimension(0:n_ene) :: energia,energia_ave !!! local variables -el - integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,& + integer :: i,ii,kkk,it,j,nacc,ijunk,indmin,indmax,& ISWEEP,Kwita,iretcode,indeold,iene,noverlap,& irep,nstart_grow,inde +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,& deix,dent,sold,scur,runtime ! @@ -1084,8 +1087,11 @@ real(kind=8),dimension(0:n_ene) :: energia,energia_ave ! !!! local variables - el - integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,& + integer :: i,j,it,ii,iproc,nacc,ISWEEP,indmax,indmin,ijunk,& Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,& frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index 87f7394..9fd36be 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -141,7 +141,9 @@ real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutine: setup_fricmat +#ifndef LBFGS real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) +#endif !----------------------------------------------------------------------------- ! ! @@ -1033,7 +1035,7 @@ ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' ! include 'COMMON.HAIRPIN' - real(kind=8),dimension(3) :: L,vcm + real(kind=8),dimension(3) :: L,vcm,boxx #ifdef VOUT real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres #endif @@ -1043,11 +1045,16 @@ !el common /gucio/ cm integer :: i,j,nharp integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3) + ! logical :: ovrtim real(kind=8) :: tt0,scalfac integer :: nres2,itime nres2=2*nres print *, "ENTER MD" + boxx(1)=boxxsize + boxx(2)=boxysize + boxx(3)=boxzsize + ! #ifdef MPI print *,"MY tmpdir",tmpdir,ilen(tmpdir) @@ -1172,7 +1179,18 @@ endif if (reset_vel .and. tbf .and. lang.eq.0 & .and. mod(itime,count_reset_vel).eq.0) then + !WARP WATER + do i=1,nres + if (molnum(i).eq.5) then + call to_box(c(1,i),c(2,i),c(3,i)) + do j=1,3 + if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j) + enddo + endif +! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i) + enddo call random_vel + write(iout,'(a,f20.2)') & "Velocities reset to random values, time",totT do i=0,2*nres @@ -2649,7 +2667,10 @@ character(len=50) :: tytul logical :: file_exist !el common /gucio/ cm - integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime + integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: etot,tt0 logical :: fail @@ -3028,6 +3049,10 @@ if (me.eq.king.or..not.out1file) write(iout,*)& 'Minimizing initial PDB structure: Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) +#ifdef LBFGS + if (me.eq.king.or..not.out1file)& + write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun +#endif else call geom_to_var(nvar,varia) if(me.eq.king.or..not.out1file) write (iout,*)& @@ -3036,9 +3061,9 @@ call var_to_geom(nvar,varia) #ifdef LBFGS if (me.eq.king.or..not.out1file)& - write(iout,*) 'LBFGS return code is ',status,' eval ',nfun + write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun if(me.eq.king.or..not.out1file)& - write(iout,*) 'LBFGS return code is ',status,' eval ',nfun + write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun #else if (me.eq.king.or..not.out1file)& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun @@ -3253,7 +3278,7 @@ !#define DEBUG #ifdef FIVEDIAG real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs - real(kind=8) :: sumx,Ek2,Ek3,aux + real(kind=8) :: sumx,Ek2,Ek3,aux,masinv #ifdef DEBUG real(kind=8) ,allocatable, dimension(:) :: rsold real (kind=8),allocatable,dimension(:,:) :: matold,inertia @@ -3290,11 +3315,15 @@ #endif do ichain=1,nchain ind=0 - if(.not.allocated(ghalf)) print *,"COCO" - if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) - ghalf=0.0d0 +! if(.not.allocated(ghalf)) print *,"COCO" +! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) +! ghalf=0.0d0 n=dimen_chain(ichain) innt=iposd_chain(ichain) + if (molnum(innt).eq.5) go to 137 + if(.not.allocated(ghalf)) print *,"COCO" + if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2)) + ghalf=0.0d0 inct=innt+n-1 #ifdef DEBUG write (iout,*) "Chain",ichain," n",n," start",innt @@ -3364,20 +3393,37 @@ enddo enddo #endif +137 continue + write(iout,*) "HERE,",n xv=0.0d0 ii=0 do i=1,n do k=1,3 ii=ii+1 + if (molnum(innt).eq.5) geigen(i)=1.0/msc(itype(innt+i-1,5),5) + sigv=dsqrt((Rb*t_bath)/geigen(i)) lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2 -!c write (iout,*) "i",i," ii",ii," geigen",geigen(i), -!c & " d_t_work_new",d_t_work_new(ii) +! write (iout,*) "i",i," ii",ii," geigen",geigen(i), & +! " d_t_work_new",d_t_work_new(ii) enddo enddo + if (molnum(innt).eq.5) then + + do k=1,3 + do i=1,n + ind=(i-1)*3+k + d_t_work(ind)=0.0d0 + masinv=1.0/msc(itype(innt+i-1,5),5) + d_t_work(ind)=d_t_work(ind)& + +masinv*d_t_work_new((i-1)*3+k) + enddo + enddo + + else do k=1,3 do i=1,n ind=(i-1)*3+k @@ -3386,10 +3432,9 @@ d_t_work(ind)=d_t_work(ind)& +Gvec(i,j)*d_t_work_new((j-1)*3+k) enddo -!c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind) -!c call flush(iout) enddo enddo + endif #ifdef DEBUG aux=0.0d0 do k=1,3 @@ -7375,7 +7420,7 @@ allocate(DM(nres2),DU1(nres2),DU2(nres2)) allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2)) !#ifdef DEBUG - allocate(Gvec(nres2,nres2)) + allocate(Gvec(1300,1300))!maximum number of protein in one chain !#endif #else write (iout,*) "Before A Allocation",nfgtasks-1 diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index 2f47cf6..6d00c0a 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -45,7 +45,7 @@ ncont=0 kkk=3 - do i=nnt+kkk,nct + do i=nnt_molec(1)+kkk,nct_molec(1) iti=iabs(itype(i,molnum(i))) if (molnum(i).lt.3) then inum=i+nres @@ -71,6 +71,7 @@ ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j) if (dist(inum,jnum).lt.rcomp) then ncont=ncont+1 + if (ncont.gt.nres*100) ncont=nres*100 icont(1,ncont)=i icont(2,ncont)=j endif @@ -172,7 +173,7 @@ ncont=0 kkk=0 ! print *,'nnt=',nnt,' nct=',nct - do i=nnt,nct-3 + do i=nnt_molec(1),nct_molec(1)-3 do k=1,3 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1)) enddo @@ -312,7 +313,7 @@ ees=0.0 evdw=0.0 ! print *, "nntt,nct",nnt,nct-2 - do 1 i=nnt,nct-2 + do 1 i=nnt_molec(1),nct_molec(1)-2 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1 xi=c(1,i) yi=c(2,i) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 516a8d7..e469077 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -1164,6 +1164,17 @@ jgrad_start(i)=i+1 jgrad_end(i)=nres enddo +! THIS SHOULD BE FIXED + itmp=0 + do i=1,4 + itmp=itmp+nres_molec(i) + enddo + call int_bounds(nres_molec(5),icatb_start,icatb_end) + icatb_start=icatb_start+itmp + icatb_end=icatb_end+itmp + + + if (lprint) then write (*,*) 'Processor:',fg_rank,' CG group',kolor,& ' absolute rank',myrank,& diff --git a/source/unres/data/energy_data.F90 b/source/unres/data/energy_data.F90 index b42fec5..8f45280 100644 --- a/source/unres/data/energy_data.F90 +++ b/source/unres/data/energy_data.F90 @@ -223,6 +223,7 @@ iphi_nucl_end,iphid_nucl_start,iphid_nucl_end,& ibondp_nucl_start,ibondp_nucl_end,ithet_nucl_start,ithet_nucl_end,& loc_start_nucl,loc_end_nucl + integer :: icatb_start,icatb_end integer,dimension(:),allocatable :: ibond_displ,ibond_count,& ithet_displ,ithet_count,iphi_displ,iphi_count,iphi1_displ,& iphi1_count,ivec_displ,ivec_count,iset_displ,iset_count,& @@ -453,7 +454,8 @@ newcontlistcatscangfi,newcontlistcatscangfj,& newcontlistcatscangfk,& newcontlistcatscangti,newcontlistcatscangtj,& - newcontlistcatscangtk,newcontlistcatscangtl + newcontlistcatscangtk,newcontlistcatscangtl,& + newcontlistcatcatnormi,newcontlistcatcatnormj @@ -468,7 +470,8 @@ g_ilist_catptran,g_ilist_catscang,g_ilist_catscangf,g_ilist_catscangt,& g_listcatscang_start,g_listcatscang_end,& g_listcatscangf_start,g_listcatscangf_end,& - g_listcatscangt_start,g_listcatscangt_end + g_listcatscangt_start,g_listcatscangt_end,& + g_listcatcatnorm_start,g_listcatcatnorm_end,g_ilist_catcatnorm ! MARTINI FORCE FIELD @@ -520,5 +523,9 @@ ! buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick !------------------------------------------------------------------------- real(kind=8),dimension(3,70000) :: ea +#ifdef LBFGS + character*9 statusbf + integer niter,nfun,ncalls +#endif ! real(kind=8) :: buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick end module energy_data diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index d81310f..572dd3e 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -274,8 +274,10 @@ integer ishield_listbuf(-1:nres_molec(1)), & shield_listbuf(maxcontsshi,-1:nres_molec(1)),k,j,i,iii,impishi,mojint,jjj + integer :: imatupdate2 ! print *,"I START ENERGY" imatupdate=100 + imatupdate2=100 ! if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list ! real(kind=8), dimension(:),allocatable:: fac_shieldbuf ! real(kind=8), dimension(:,:,:),allocatable:: & @@ -409,7 +411,7 @@ time_Bcastw=time_Bcastw+MPI_Wtime()-time00 ! call chainbuild_cart endif - print *,"itime_mat",itime_mat,imatupdate +! print *,"itime_mat",itime_mat,imatupdate if (nfgtasks.gt.1) then call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) endif @@ -424,9 +426,18 @@ if (mod(itime_mat,imatupdate).eq.0) then ! print *,'Processor',myrank,' calling etotal ipot=',ipot call make_cat_pep_list +! call make_cat_cat_list endif endif endif + if (nres_molec(5).gt.0) then + if (mod(itime_mat,imatupdate2).eq.0) then +! print *, "before cat cat" +! print *,'Processor',myrank,' calling etotal ipot=',ipot +! call make_cat_pep_list + call make_cat_cat_list + endif + endif ! write (iout,*) "after make_pp_inter_list" ! print *,'Processor',myrank,' calling etotal ipot=',ipot @@ -893,7 +904,7 @@ Eafmforce=0.0d0 endif endif - print *,"before tubemode",tubemode +! print *,"before tubemode",tubemode if (tubemode.eq.1) then call calctube(etube) else if (tubemode.eq.2) then @@ -903,7 +914,7 @@ else etube=0.0d0 endif - print *, "TU JEST PRZED EHPB" +! print *, "TU JEST PRZED EHPB" call edis(ehpb) !-------------------------------------------------------- @@ -949,13 +960,13 @@ else ecation_protang=0.0d0 endif - if (nfgtasks.gt.1) then - if (fg_rank.eq.0) then +! if (nfgtasks.gt.1) then +! if (fg_rank.eq.0) then if (nres_molec(5).gt.1) call ecatcat(ecationcation) - endif - else - if (nres_molec(5).gt.1) call ecatcat(ecationcation) - endif +! endif +! else +! if (nres_molec(5).gt.1) call ecatcat(ecationcation) +! endif if (oldion.gt.0) then if (g_ilist_catpnorm.gt.0) call ecat_prot(ecation_prot) else @@ -13420,8 +13431,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! include 'COMMON.VAR' ! include 'COMMON.CONTACTS' use comm_srutu +!#ifdef LBFGS +! use minimm, only: funcgrad +!#endif !el integer :: icall !el common /srutu/ icall +! real(kind=8) :: funcgrad real(kind=8),dimension(6) :: ggg real(kind=8),dimension(3) :: cc,xx,ddc,ddx real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) @@ -13431,7 +13446,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! real(kind=8) :: urparm(1) !EL external fdum integer :: nf,i,j,k - real(kind=8) :: aincr,etot,etot1 + real(kind=8) :: aincr,etot,etot1,ff icg=1 nf=0 nfl=0 @@ -13443,8 +13458,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! call geom_to_var(nvar,x) call etotal(energia) etot=energia(0) +#ifdef LBFGS + ff=funcgrad(x,g) +#else !el call enerprint(energia) call gradient(nvar,x,nf,g,uiparm,urparm,fdum) +#endif icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) @@ -13997,8 +14016,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! include 'COMMON.VAR' ! include 'COMMON.GEO' use comm_srutu +!#ifdef LBFGS +! use minimm, only : funcgrad +!#endif !el integer :: icall !el common /srutu/ icall +! real(kind=8) :: funcgrad real(kind=8),dimension(6*nres) :: x,gana,gg !(maxvar) (maxvar=6*maxres) integer :: uiparm(1) real(kind=8) :: urparm(1) @@ -14006,7 +14029,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! character(len=6) :: key !EL external fdum integer :: i,ii,nf - real(kind=8) :: xi,aincr,etot,etot1,etot2 + real(kind=8) :: xi,aincr,etot,etot1,etot2,ff call zerograd aincr=1.0D-7 print '(a)','Calling CHECK_INT.' @@ -14032,9 +14055,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #endif nf=1 nfl=3 +#ifdef LBFGS + ff=funcgrad(x,gana) +#else + !d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) !d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp +#endif icall=1 do i=1,nvar xi=x(i) @@ -17713,6 +17741,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !----------------------------------------------------------------------------- ! gradient_p.F !----------------------------------------------------------------------------- +#ifndef LBFGS subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm) use io_base, only:intout,briefout @@ -17818,6 +17847,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !d write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n) return end subroutine gradient +#endif !----------------------------------------------------------------------------- subroutine func(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F @@ -21972,6 +22002,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! allocate(newcontlistcatscnormj(300*nres)) allocate(newcontlistcatpnormi(300*nres)) allocate(newcontlistcatpnormj(300*nres)) + allocate(newcontlistcatcatnormi(900*nres)) + allocate(newcontlistcatcatnormj(900*nres)) allocate(newcontlistcatscangi(300*nres)) allocate(newcontlistcatscangj(300*nres)) @@ -23608,7 +23640,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #endif subroutine ecatcat(ecationcation) use MD_data, only: t_bath - integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff + integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff,& + ii real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& r7,r4,ecationcation,k0,rcal,aa,bb,sslipi,ssgradlipi,sslipj,ssgradlipj real(kind=8) :: xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & @@ -23626,12 +23659,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! k0 = 332.0*(2.0*2.0)/80.0 itmp=0 - do i=1,4 - itmp=itmp+nres_molec(i) - enddo -! write(iout,*) "itmp",itmp - do i=itmp+1,itmp+nres_molec(5)-1 - +! do i=1,4 +! itmp=itmp+nres_molec(i) +! enddo +! write(iout,*) "itmp",g_listcatcatnorm_start, g_listcatcatnorm_end +! do i=itmp+1,itmp+nres_molec(5)-1 + do ii=g_listcatcatnorm_start, g_listcatcatnorm_end + i=newcontlistcatcatnormi(ii) + j=newcontlistcatcatnormj(ii) + xi=c(1,i) yi=c(2,i) zi=c(3,i) @@ -23639,7 +23675,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! itypi=itype(i,5) call to_box(xi,yi,zi) call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) - do j=i+1,itmp+nres_molec(5) +! do j=i+1,itmp+nres_molec(5) itypj=itype(j,5) ! print *,i,j,itypi,itypj k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0 @@ -23659,7 +23695,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! rcal =xj**2+yj**2+zj**2 ract=sqrt(rcal) if ((itypi.gt.1).or.(itypj.gt.1)) then - + if (sss2min2.eq.0.0d0) cycle + sss2min2=sscale2(ract,12.0d0,1.0d0) + sss2mingrad2=sscagrad2(ract,12.0d0,1.0d0) ! rcat0=3.472 ! epscalc=0.05 ! r06 = rcat0**6 @@ -23680,13 +23718,13 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! enddo do k=1,3 gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k) - gradcatcat(k,i)=gradcatcat(k,i)-gg(k) - gradcatcat(k,j)=gradcatcat(k,j)+gg(k) + gradcatcat(k,i)=gradcatcat(k,i)-(gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2) + gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2 enddo if (energy_dec) write (iout,*) "ecatcat",i,j,Evan1cat,Evan2cat,Eeleccat,& r012,rcal**6,ichargecat(itypi)*ichargecat(itypj) ! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj - ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat + ecationcation=ecationcation+(Evan1cat+Evan2cat+Eeleccat)*sss2min2 else !this is water part and other non standard molecules sss2min2=sscale2(ract,10.0d0,1.0d0)! cutoff for water interaction is 15A @@ -23715,10 +23753,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! gradcatcat(k,i)=gradcatcat(k,i)-gg(k)*sss2min2-sss2mingrad2*ewater*r(k)/ract gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+sss2mingrad2*ewater*r(k)/ract enddo - if (energy_dec) write(iout,'(2f10.7,f15.7,2i5)') rdiff,ract,ecationcation,i,j + if (energy_dec) write(iout,'(2f8.2,f10.2,2i5)') rdiff,ract,ecationcation,i,j endif ! end water enddo - enddo +! enddo return end subroutine ecatcat !--------------------------------------------------------------------------- @@ -29718,7 +29756,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! real(kind=8) :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi real(kind=8) :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj real(kind=8) :: xja,yja,zja - integer:: contlistcatpnormi(250*nres),contlistcatpnormj(250*nres) + integer:: contlistcatpnormi(300*nres),contlistcatpnormj(300*nres) integer:: contlistcatscnormi(250*nres),contlistcatscnormj(250*nres) integer:: contlistcatptrani(250*nres),contlistcatptranj(250*nres) integer:: contlistcatsctrani(250*nres),contlistcatsctranj(250*nres) @@ -30162,6 +30200,134 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! return end subroutine make_cat_pep_list + subroutine make_cat_cat_list + include 'mpif.h' + real(kind=8) :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp + real(kind=8) :: xmedj,ymedj,zmedj,sslipi,ssgradlipi,faclipij2,sslipj,ssgradlipj + real(kind=8) :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi + real(kind=8) :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj + real(kind=8) :: xja,yja,zja + integer,dimension(:),allocatable:: contlistcatpnormi,contlistcatpnormj +! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_catscnorm,& + ilist_catsctran,ilist_catpnorm,ilist_catptran,itmp,ilist_catscang,& + ilist_catscangf,ilist_catscangt,k + integer displ(0:nprocs),i_ilist_catscnorm(0:nprocs),ierr,& + i_ilist_catpnorm(0:nprocs),i_ilist_catsctran(0:nprocs),& + i_ilist_catptran(0:nprocs),i_ilist_catscang(0:nprocs),& + i_ilist_catscangf(0:nprocs),i_ilist_catscangt(0:nprocs) + write(iout,*),"START make_catcat" + ilist_catpnorm=0 + ilist_catscnorm=0 + ilist_catptran=0 + ilist_catsctran=0 + ilist_catscang=0 + + if (.not.allocated(contlistcatpnormi)) then + allocate(contlistcatpnormi(900*nres)) + allocate(contlistcatpnormj(900*nres)) + endif + r_buff_list=3.0 + itmp=0 + do i=1,4 + itmp=itmp+nres_molec(i) + enddo +! go to 17 +! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization + do i=icatb_start,icatb_end + xi=c(1,i) + yi=c(2,i) + zi=c(3,i) + call to_box(xi,yi,zi) + dxi=dc_norm(1,i) + dyi=dc_norm(2,i) + dzi=dc_norm(3,i) +! dsci_inv=vbld_inv(i+nres) + do j=i+1,itmp+nres_molec(5) + dxj=dc(1,j) + dyj=dc(2,j) + dzj=dc(3,j) + dx_normj=dc_norm(1,j) + dy_normj=dc_norm(2,j) + dz_normj=dc_norm(3,j) + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + call to_box(xj,yj,zj) +! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) +! faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 + xja=boxshift(xj-xi,boxxsize) + yja=boxshift(yj-yi,boxysize) + zja=boxshift(zj-zi,boxzsize) + dist_init=xja**2+yja**2+zja**2 + if (sqrt(dist_init).le.(10.0+r_buff_list)) then +! Here the list is created +! if (i.eq.2) then +! print *,i,j,dist_init,ilist_catpnorm +! endif + ilist_catpnorm=ilist_catpnorm+1 + +! this can be substituted by cantor and anti-cantor + contlistcatpnormi(ilist_catpnorm)=i + contlistcatpnormj(ilist_catpnorm)=j + endif +! enddo + enddo + enddo +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_catsctran,ilist_catptran,& + ilist_catscnorm,ilist_catpnorm,ilist_catscang + + do i=1,ilist_catpnorm + write (iout,*) i,contlistcatpnormi(i) + enddo + + +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_catpnorm,g_ilist_catcatnorm,1,& + MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +! write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_catpnorm,1,MPI_INTEGER,& + i_ilist_catpnorm,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_catpnorm(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistcatpnormi,ilist_catpnorm,MPI_INTEGER,& + newcontlistcatcatnormi,i_ilist_catpnorm,displ,MPI_INTEGER,& + king,FG_COMM,IERR) + call MPI_Gatherv(contlistcatpnormj,ilist_catpnorm,MPI_INTEGER,& + newcontlistcatcatnormj,i_ilist_catpnorm,displ,MPI_INTEGER,& + king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_catcatnorm,1,MPI_INT,king,FG_COMM,IERR) +! write(iout,*) "before bcast",g_ilist_sc +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + call MPI_Bcast(newcontlistcatcatnormi,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR) + call MPI_Bcast(newcontlistcatcatnormj,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR) + + + else + g_ilist_catcatnorm=ilist_catpnorm + do i=1,ilist_catpnorm + newcontlistcatcatnormi(i)=contlistcatpnormi(i) + newcontlistcatcatnormj(i)=contlistcatpnormj(i) + enddo + endif + call int_bounds(g_ilist_catcatnorm,g_listcatcatnorm_start,g_listcatcatnorm_end) + +#ifdef DEBUG + write (iout,*) "after MPIREDUCE",g_ilist_catcatnorm + + do i=1,g_ilist_catcatnorm + write (iout,*) i,newcontlistcatcatnormi(i),newcontlistcatcatnormj(i) + enddo +#endif + write(iout,*),"END make_catcat" + return + end subroutine make_cat_cat_list !----------------------------------------------------------------------------- @@ -31832,5 +31998,255 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! return end function sq +#ifdef LBFGS + double precision function funcgrad(x,g) + use MD_data, only: totT,usampl + implicit none + double precision energia(0:n_ene) + double precision x(nvar),g(nvar) + integer i + call var_to_geom(nvar,x) + call zerograd + call chainbuild + call etotal(energia(0)) + call sum_gradient + funcgrad=energia(0) + call cart2intgrad(nvar,g) + 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 + g(i)=g(i)+gloc(i,icg) + enddo + return + end function funcgrad + subroutine cart2intgrad(n,g) + integer n + double precision g(n) + double precision drt(3,3,nres),rdt(3,3,nres),dp(3,3),& + temp(3,3),prordt(3,3,nres),prodrt(3,3,nres) + double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp + double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,& + cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl + double precision fromto(3,3),aux(6) + integer i,ii,j,jjj,k,l,m,indi,ind,ind1 + logical sideonly + sideonly=.false. + g=0.0d0 + if (sideonly) goto 10 + do i=1,nres-2 + rdt(1,1,i)=-rt(1,2,i) + rdt(1,2,i)= rt(1,1,i) + rdt(1,3,i)= 0.0d0 + rdt(2,1,i)=-rt(2,2,i) + rdt(2,2,i)= rt(2,1,i) + rdt(2,3,i)= 0.0d0 + rdt(3,1,i)=-rt(3,2,i) + rdt(3,2,i)= rt(3,1,i) + rdt(3,3,i)= 0.0d0 + enddo + do i=2,nres-2 + drt(1,1,i)= 0.0d0 + drt(1,2,i)= 0.0d0 + drt(1,3,i)= 0.0d0 + drt(2,1,i)= rt(3,1,i) + drt(2,2,i)= rt(3,2,i) + drt(2,3,i)= rt(3,3,i) + drt(3,1,i)=-rt(2,1,i) + drt(3,2,i)=-rt(2,2,i) + drt(3,3,i)=-rt(2,3,i) + enddo + ind1=0 + do i=1,nres-2 + ind1=ind1+1 + if (n.gt.nphi) then + + do j=1,3 + do k=1,2 + dpjk=0.0D0 + do l=1,3 + dpjk=dpjk+prod(j,l,i)*rdt(l,k,i) + enddo + dp(j,k)=dpjk + prordt(j,k,i)=dp(j,k) + enddo + dp(j,3)=0.0D0 + g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg) + enddo + xx1(1)=-0.5D0*xloc(2,i+1) + xx1(2)= 0.5D0*xloc(1,i+1) + do j=1,3 + xj=0.0D0 + do k=1,2 + xj=xj+r(j,k,i)*xx1(k) + enddo + xx(j)=xj + enddo + do j=1,3 + rj=0.0D0 + do k=1,3 + rj=rj+prod(j,k,i)*xx(k) + enddo + g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg) + enddo + if (i.lt.nres-2) then + do j=1,3 + dxoiij=0.0D0 + do k=1,3 + dxoiij=dxoiij+dp(j,k)*xrot(k,i+2) + enddo + g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg) + enddo + endif + + endif + + + if (i.gt.1) then + do j=1,3 + do k=1,3 + dpjk=0.0 + do l=2,3 + dpjk=dpjk+prod(j,l,i)*drt(l,k,i) + enddo + dp(j,k)=dpjk + prodrt(j,k,i)=dp(j,k) + enddo + g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg) + enddo + endif + xx(1)= 0.0D0 + xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i) + xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i) + if (i.gt.1) then + do j=1,3 + rj=0.0D0 + do k=2,3 + rj=rj+prod(j,k,i)*xx(k) + enddo + g(i-1)=g(i-1)-rj*gradx(j,i+1,icg) + enddo + endif + if (i.gt.1) then + do j=1,3 + dxoiij=0.0D0 + do k=1,3 + dxoiij=dxoiij+dp(j,k)*xrot(k,i+2) + enddo + g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg) + enddo + endif + do j=i+1,nres-2 + ind1=ind1+1 + call build_fromto(i+1,j+1,fromto) + do k=1,3 + do l=1,3 + tempkl=0.0D0 + do m=1,2 + tempkl=tempkl+prordt(k,m,i)*fromto(m,l) + enddo + temp(k,l)=tempkl + enddo + enddo + if (n.gt.nphi) then + do k=1,3 + g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg) + enddo + do k=1,3 + dxoijk=0.0D0 + do l=1,3 + dxoijk=dxoijk+temp(k,l)*xrot(l,j+2) + enddo + g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg) + enddo + endif + do k=1,3 + do l=1,3 + tempkl=0.0D0 + do m=1,3 + tempkl=tempkl+prodrt(k,m,i)*fromto(m,l) + enddo + temp(k,l)=tempkl + enddo + enddo + if (i.gt.1) then + do k=1,3 + g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg) + enddo + do k=1,3 + dxoijk=0.0D0 + do l=1,3 + dxoijk=dxoijk+temp(k,l)*xrot(l,j+2) + enddo + g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg) + enddo + endif + enddo + enddo + + if (nvar.le.nphi+ntheta) return + + 10 continue + do i=2,nres-1 + if (iabs(itype(i,1)).eq.10 .or. itype(i,1).eq.ntyp1& !) cycle + .or. mask_side(i).eq.0 ) cycle + ii=ialph(i,1) + dsci=vbld(i+nres) +#ifdef OSF + alphi=alph(i) + omegi=omeg(i) + if(alphi.ne.alphi) alphi=100.0 + if(omegi.ne.omegi) omegi=-100.0 +#else + alphi=alph(i) + omegi=omeg(i) +#endif + cosalphi=dcos(alphi) + sinalphi=dsin(alphi) + cosomegi=dcos(omegi) + sinomegi=dsin(omegi) + temp(1,1)=-dsci*sinalphi + temp(2,1)= dsci*cosalphi*cosomegi + temp(3,1)=-dsci*cosalphi*sinomegi + temp(1,2)=0.0D0 + temp(2,2)=-dsci*sinalphi*sinomegi + temp(3,2)=-dsci*sinalphi*cosomegi + theta2=pi-0.5D0*theta(i+1) + cost2=dcos(theta2) + sint2=dsin(theta2) + jjj=0 + do j=1,2 + xp=temp(1,j) + yp=temp(2,j) + xxp= xp*cost2+yp*sint2 + yyp=-xp*sint2+yp*cost2 + zzp=temp(3,j) + xx(1)=xxp + xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1) + xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1) + do k=1,3 + dj=0.0D0 + do l=1,3 + dj=dj+prod(k,l,i-1)*xx(l) + enddo + aux(jjj+k)=dj + enddo + jjj=jjj+3 + enddo + do k=1,3 + g(ii)=g(ii)+aux(k)*gradx(k,i,icg) + g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg) + enddo + enddo + return + end subroutine cart2intgrad + + +#endif !-------------------------------------------------------------------------- end module energy diff --git a/source/unres/io.F90 b/source/unres/io.F90 index 4750e57..3f3c127 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -278,7 +278,7 @@ enddo enddo - call contact(.true.,ncont_ref,icont_ref,co) + call contact(.false.,ncont_ref,icont_ref,co) ! do k=1,nres ! write(icsa_native_int,200) (crefjlee(i,k),i=1,3) @@ -1464,7 +1464,7 @@ cref(j,i,kkk)=c(j,i) enddo enddo - call contact(.true.,ncont_ref,icont_ref,co) + call contact(.false.,ncont_ref,icont_ref,co) endif ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup ! call flush(iout) @@ -1477,6 +1477,7 @@ if (pdbref) then if(me.eq.king.or..not.out1file) & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup + if (.false.) then do i=1,ncont_ref do j=1,2 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup @@ -1487,6 +1488,7 @@ restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i) enddo endif + endif if (constr_homology.gt.0) then ! write (iout,*) "Calling read_constr_homology" ! call flush(iout) @@ -2839,7 +2841,7 @@ write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain) enddo return - end + end subroutine !c--------------------------------------------------------------------- integer function tperm(i,iperm,tabpermchain) ! implicit none diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90 index 1d7c037..f6058d7 100644 --- a/source/unres/io_base.F90 +++ b/source/unres/io_base.F90 @@ -155,9 +155,9 @@ !----------------------------------------------------------------------------- subroutine read_x(kanal,*) - use geometry_data - use energy_data - use geometry, only:int_from_cart1 + use geometry_data + use energy_data + use geometry, only:int_from_cart1 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' @@ -1018,7 +1018,6 @@ use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize use energy_data -! use energy, only: to_box ! use control use compare_data use MD_data @@ -1142,7 +1141,7 @@ ires=0 write(iout,*) "TUTUT" do i=nnt,nct - write(iout,*), "coord",c(1,i),c(2,i),c(3,i) +! write(iout,*), "coord",c(1,i),c(2,i),c(3,i) iti=itype(i,molnum(i)) print *,i,molnum(i) if (molnum(i+1).eq.0) then @@ -1183,7 +1182,20 @@ ! if (zi.lt.0.0d0) zi=zi+boxzsize write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,xi,yi,zi,vtot(i) - else + elseif (molnum(i).eq.5) then + xi=c(1,i) + yi=c(2,i) + zi=c(3,i) + xi=dmod(xi,boxxsize) + if (xi.lt.0.0d0) xi=xi+boxxsize + yi=dmod(yi,boxysize) + if (yi.lt.0.0d0) yi=yi+boxysize + zi=dmod(zi,boxzsize) + if (zi.lt.0.0d0) zi=zi+boxzsize + write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& + ires,xi,yi,zi,vtot(i) + else + write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) endif @@ -1204,6 +1216,7 @@ enddo write (iunit,'(a)') 'TER' do i=nnt,nct-1 + if (molnum(i).eq.5) cycle if (itype(i,1).eq.ntyp1) cycle if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then write (iunit,30) ica(i),ica(i+1) @@ -1560,7 +1573,6 @@ !----------------------------------------------------------------------------- integer function rescode(iseq,nam,itype,molecule) -! use io_base, only: ucase ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.NAMES' diff --git a/source/unres/map.F90 b/source/unres/map.F90 index b91d43e..37aa22e 100644 --- a/source/unres/map.F90 +++ b/source/unres/map.F90 @@ -91,12 +91,15 @@ real(kind=8) :: energia(0:n_ene) character(len=5) :: angid(4)=reshape((/'PHI ','THETA','ALPHA',& 'OMEGA'/),shape(angid)) - real(kind=8) :: ang_list(10) + real(kind=8) :: ang_list(10),ff !el real(kind=8),dimension(6*nres) :: g,x !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(:),allocatable :: g,x !(maxvar) (maxvar=6*maxres) integer :: nn(10) !el local variables - integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode,nfun + integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: etot,gnorm!,fdum integer,dimension(1) :: uiparm real(kind=8),dimension(1) :: urparm @@ -157,7 +160,14 @@ etot = energia(0) nf=1 nfl=3 - call gradient(nvar,x,nf,g,uiparm,urparm,fdum) +#ifdef LBFGS + ff=funcgrad(x,g) +#else + +!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) + call gradient(nvar,x,nf,g,uiparm,urparm,fdum) +!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp +#endif gnorm=0.0d0 do k=1,nvar gnorm=gnorm+g(k)**2 diff --git a/source/unres/minim.F90 b/source/unres/minim.F90 index 8a4020f..f15d802 100644 --- a/source/unres/minim.F90 +++ b/source/unres/minim.F90 @@ -12,6 +12,9 @@ ! use csa_data ! use energy implicit none +!#ifdef LBFGS +! double precision funcgrad_dc +!#endif !----------------------------------------------------------------------------- ! ! @@ -25,7 +28,7 @@ ! *** assess candidate step (***sol version 2.3) *** ! integer :: liv, l - integer(kind=8) ::lv + integer ::lv integer :: iv(liv) real(kind=8) :: v(lv) ! @@ -545,7 +548,7 @@ ! *** alg = 2 means general unconstrained optimization constants. ! integer :: liv, l - integer(kind=8)::lv + integer::lv integer :: iv(liv) integer :: alg real(kind=8) :: v(lv) @@ -705,7 +708,7 @@ ! *** parameter declarations *** ! integer :: liv, p - integer(kind=8)::lv + integer::lv integer :: iv(liv) real(kind=8) :: d(p), g(p), v(lv), x(p) ! @@ -1041,7 +1044,7 @@ ! *** alg = 1 for regression, alg = 2 for general unconstrained opt. ! integer :: alg, liv, n - integer(kind=8) :: lv + integer :: lv integer :: iv(liv) real(kind=8) :: d(n), v(lv) ! @@ -1430,7 +1433,7 @@ ! *** alg = 2 means general unconstrained optimization constants. ! integer :: alg, l - integer(kind=8)::lv + integer::lv real(kind=8) :: v(lv) !/+ !el real(kind=8) :: dmax1 @@ -1612,7 +1615,7 @@ ! *** (analytic) gradient and hessian provided by the caller. *** ! integer :: liv, n - integer(kind=8):: lv + integer:: lv integer:: iv(liv) integer :: uiparm(1) real(kind=8) :: d(n), x(n), v(lv), urparm(1) @@ -1760,7 +1763,7 @@ ! *** parameter declarations *** ! integer :: lh, liv, n - integer(kind=8)::lv + integer::lv integer :: iv(liv) real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n) ! @@ -2210,7 +2213,7 @@ ! *** parameter declarations *** ! integer :: liv, n - integer(kind=8) ::lv + integer ::lv integer :: iv(liv) real(kind=8) :: d(n), hdiag(n), v(lv) ! @@ -3186,8 +3189,16 @@ ! minimize_p.F !----------------------------------------------------------------------------- subroutine minimize(etot,x,iretcode,nfun) - +#ifdef LBFGS + use minima + use inform + use output + use iounit + use scales + use energy, only: funcgrad,etotal,enerprint +#else use energy, only: func,gradient,fdum!,etotal,enerprint +#endif use comm_srutu ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -3214,19 +3225,54 @@ real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv) !el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) - real(kind=8) :: energia(0:n_ene) + real(kind=8) :: energia(0:n_ene),grdmin ! external func,gradient,fdum ! external func_restr,grad_restr logical :: not_done,change,reduce !el common /przechowalnia/ v !el local variables integer :: iretcode,nfun,nvar_restr,idum(1),j - integer (kind=8):: lv + integer :: lv real(kind=8) :: etot,rdum(1) !,fdum - +#ifdef LBFGS + external optsave +#endif lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) - +#ifndef LBFGS if (.not.allocated(v)) allocate(v(1:lv)) +#endif +#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 + 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 @@ -3328,12 +3374,13 @@ ! write (*,*) 'Processor',MyID,' leaves MINIMIZE.' ! ENDDO ! NOT_DONE - +#endif return end subroutine minimize !----------------------------------------------------------------------------- ! gradient_p.F !----------------------------------------------------------------------------- +#ifndef LBFGS subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm) use energy, only: cartder,zerograd,etotal,sum_gradient @@ -3489,6 +3536,8 @@ !d enddo return end subroutine grad_restr +#endif +#ifndef LBFGS !----------------------------------------------------------------------------- subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F @@ -3528,6 +3577,7 @@ ! endif return end subroutine func_restr +#endif !----------------------------------------------------------------------------- ! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy !----------------------------------------------------------------------------- @@ -3585,9 +3635,16 @@ !el subroutine xx2x(x,xx) in module math !----------------------------------------------------------------------------- subroutine minim_dc(etot,iretcode,nfun) +#ifdef LBFGS + use minima + use inform + use output + use iounit + use scales +#endif use MPI_data - use energy, only: fdum,check_ecartint + use energy, only: fdum,check_ecartint,etotal,enerprint ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI @@ -3602,16 +3659,26 @@ ! include 'COMMON.MINIM' ! include 'COMMON.CHAIN' integer :: iretcode,nfun,k,i,j,idum(1) - integer(kind=8) :: lv + integer :: lv,nvarx,nf integer,dimension(liv) :: iv real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv) - real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres) + real(kind=8),dimension(:),allocatable :: x,d,xx !(maxvar) (maxvar=6*maxres) !el common /przechowalnia/ v - real(kind=8) :: energia(0:n_ene) + real(kind=8) :: energia(0:n_ene),grdmin ! external func_dc,grad_dc ,fdum logical :: not_done,change,reduce real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum +#ifdef LBFGS + external optsave + maxiter=maxmin + coordtype='CARTESIAN' + grdmin=tolf + jout=iout + jprint=print_min_stat + iwrite=0 + write(iout,*) "in minim_dc LBFGS" +#else lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) print *,"lv value",lv @@ -3661,7 +3728,12 @@ do i=1,6*nres d(i)=1.0D-1 enddo - +#endif + if (.not.allocated(x)) then + allocate(x(6*nres)) + allocate(xx(6*nres)) + allocate(d(6*nres)) + endif k=0 do i=1,nres-1 do j=1,3 @@ -3677,9 +3749,39 @@ 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 + write (iout,*) "minim_dc Calling lbfgs" + call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave) + deallocate(scale) +#else ! call check_ecartint call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum) ! call check_ecartint +#endif k=0 do i=1,nres-1 do j=1,3 @@ -3719,11 +3821,68 @@ ! (c(j,i+nres),j=1,3) ! enddo !el---------------------------- +#ifndef LBFGS etot=v(10) iretcode=iv(1) nfun=iv(6) +#endif return end subroutine minim_dc + + +#ifdef LBFGS + real(kind=8) function funcgrad_dc(x,g) + use energy, only: etotal, zerograd,cartgrad +! implicit real*8 (a-h,o-z) +#ifdef MPI + include 'mpif.h' +#endif + integer k,nf,i,j + real(kind=8),dimension(nres*6):: x,g + 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) + 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 function +#else !----------------------------------------------------------------------------- subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm) @@ -3858,6 +4017,7 @@ return end subroutine grad_dc +#endif !----------------------------------------------------------------------------- ! minim_mcmf.F !----------------------------------------------------------------------------- @@ -3866,6 +4026,7 @@ use MPI_data use csa_data +#ifndef LBFGS use energy, only: func,gradient,fdum ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -3889,7 +4050,7 @@ integer,dimension(6) :: indx integer,dimension(liv) :: iv integer :: idum(1),nf ! - integer(kind=8) :: lv + integer :: lv real(kind=8) :: rdum(1) real(kind=8) :: przes(3),obrot(3,3),eee logical :: non_conv @@ -3995,6 +4156,7 @@ call mpi_send(ene0,1,mpi_double_precision,king,idreal,& CG_COMM,ierr) go to 10 +#endif return end subroutine minim_mcmf #endif @@ -4595,8 +4757,9 @@ end subroutine sc_minimize !----------------------------------------------------------------------------- subroutine minimize_sc1(etot,x,iretcode,nfun) - +#ifndef LBFGS use energy, only: func,gradient,fdum,etotal,enerprint +#endif use comm_srutu ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -4618,11 +4781,11 @@ ! external func_restr1,grad_restr1 logical :: not_done,change,reduce !el common /przechowalnia/ v - integer(kind=8):: lv + integer:: lv integer :: nvar_restr,i,j integer :: idum(1) real(kind=8) :: rdum(1),etot !,fdum - +#ifndef LBFGS lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres) if (.not. allocated(v)) allocate(v(1:lv)) @@ -4697,7 +4860,7 @@ etot=v(10) iretcode=iv(1) nfun=iv(6) - +#endif return end subroutine minimize_sc1 !----------------------------------------------------------------------------- @@ -4762,6 +4925,7 @@ return end subroutine func_restr1 !----------------------------------------------------------------------------- + subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm) use energy, only: cartder,zerograd,esc,sc_grad @@ -5061,7 +5225,7 @@ ! ! use control integer :: n, liv - integer(kind=8)::lv + integer::lv integer:: iv(liv) integer :: uiparm(1) real(kind=8) :: d(n), x(n), v(lv), urparm(1) @@ -5491,7 +5655,7 @@ ! integer :: iv1, nf real(kind=8) :: f - integer(kind=8)::g1 + integer::g1 ! ! *** subscripts for iv *** ! @@ -5565,7 +5729,7 @@ ! *** parameter declarations *** ! integer :: liv, n - integer(kind=8) :: lv + integer :: lv integer :: iv(liv) real(kind=8) :: d(n), fx, g(n), v(lv), x(n) ! @@ -6020,7 +6184,7 @@ ! *** parameter declarations *** ! integer :: n - integer(kind=8) :: lv + integer :: lv real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv) ! ! *** purpose *** @@ -6539,5 +6703,501 @@ ! *** last card of wzbfgs follows *** end subroutine wzbfgs !----------------------------------------------------------------------------- + +! +! +! ################################################### +! ## COPYRIGHT (C) 1999 by Jay William Ponder ## +! ## All Rights Reserved ## +! ################################################### +! +! ############################################################## +! ## ## +! ## subroutine lbfgs -- limited memory BFGS optimization ## +! ## ## +! ############################################################## +! +! +! "lbfgs" is a limited memory BFGS quasi-newton nonlinear +! optimization routine +! +! literature references: +! +! J. Nocedal, "Updating Quasi-Newton Matrices with Limited +! Storage", Mathematics of Computation, 35, 773-782 (1980) +! +! D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method +! for Large Scale Optimization", Mathematical Programming, +! 45, 503-528 (1989) +! +! J. Nocedal and S. J. Wright, "Numerical Optimization", +! Springer-Verlag, New York, 1999, Section 9.1 +! +! variables and parameters: +! +! nvar number of parameters in the objective function +! x0 contains starting point upon input, upon return +! contains the best point found +! minimum during optimization contains best current function +! value; returns final best function value +! grdmin normal exit if rms gradient gets below this value +! ncalls total number of function/gradient evaluations +! +! required external routines: +! +! fgvalue function to evaluate function and gradient values +! optsave subroutine to write out info about current status +! +! + subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) + use inform + use iounit + use keys + use linmin + use math2 + use minima + use output + use scales + use MD_data, only: itime_mat +#ifdef LBFGS + use energy_data,only:statusbf,niter,ncalls +#endif + implicit none + integer i,j,k,m + integer nvar,next + integer msav,muse +#ifndef LBFGS + integer niter,ncalls + character*9 statusbf +#endif + integer nerr,maxerr + real*8 f,f_old,fgvalue + real*8 f_move,x_move + real*8 g_norm,g_rms + real*8 minimum,grdmin + real*8 angle,rms,beta + real*8 ys,yy,gamma + real*8 x0(*) + real*8, allocatable :: rho(:) + real*8, allocatable :: alpha(:) + real*8, allocatable :: x_old(:) + real*8, allocatable :: g(:) + real*8, allocatable :: g_old(:) + real*8, allocatable :: p(:) + real*8, allocatable :: q(:) + real*8, allocatable :: r(:) + real*8, allocatable :: h0(:) + real*8, allocatable :: s(:,:) + real*8, allocatable :: y(:,:) + logical done + character*9 blank + character*20 keyword + character*240 record + character*240 string + external fgvalue,optsave +! common /lbfgstat/ status,niter,ncalls +!c +!c +!c initialize some values to be used below +!c + ncalls = 0 + rms = sqrt(dble(nvar)) + if (coordtype .eq. 'CARTESIAN') then + rms = rms / sqrt(3.0d0) + else if (coordtype .eq. 'RIGIDBODY') then + rms = rms / sqrt(6.0d0) + end if + blank = ' ' + done = .false. + nerr = 0 + maxerr = 2 +!c +!c perform dynamic allocation of some global arrays +!c + if (.not. allocated(scale)) allocate (scale(nvar)) +!c +!c set default values for variable scale factors +!c + if (.not. set_scale) then + do i = 1, nvar + if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0 + end do + end if +!c +!c set default parameters for the optimization +!c + msav = min(nvar,20) + if (fctmin .eq. 0.0d0) fctmin = -100000000.0d0 + if (maxiter .eq. 0) maxiter = 1000000 + if (nextiter .eq. 0) nextiter = 1 + if (jprint .lt. 0) jprint = 1 + if (iwrite .lt. 0) iwrite = 1 + write(iout,*) "maxiter",maxiter +!c +!c set default parameters for the line search +!c + if (stpmax .eq. 0.0d0) stpmax = 5.0d0 + stpmin = 1.0d-16 + cappa = 0.9d0 + slpmax = 10000.0d0 + angmax = 180.0d0 + intmax = 5 +!c +!c search the keywords for optimization parameters +!c +#ifdef LBFGSREAD + do i = 1, nkey + next = 1 + record = keyline(i) + call gettext (record,keyword,next) + call upcase (keyword) + string = record(next:240) + if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then + read (string,*,err=10,end=10) msav + msav = max(0,min(msav,nvar)) + else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then + msav = 0 + else if (keyword(1:7) .eq. 'FCTMIN ') then + read (string,*,err=10,end=10) fctmin + else if (keyword(1:8) .eq. 'MAXITER ') then + read (string,*,err=10,end=10) maxiter + else if (keyword(1:8) .eq. 'STEPMAX ') then + read (string,*,err=10,end=10) stpmax + else if (keyword(1:8) .eq. 'STEPMIN ') then + read (string,*,err=10,end=10) stpmin + else if (keyword(1:6) .eq. 'CAPPA ') then + read (string,*,err=10,end=10) cappa + else if (keyword(1:9) .eq. 'SLOPEMAX ') then + read (string,*,err=10,end=10) slpmax + else if (keyword(1:7) .eq. 'ANGMAX ') then + read (string,*,err=10,end=10) angmax + else if (keyword(1:7) .eq. 'INTMAX ') then + read (string,*,err=10,end=10) intmax + end if + 10 continue + end do +#endif +!c +!c print header information about the optimization method +!c + if (jprint .gt. 0) then + if (msav .eq. 0) then + write (jout,20) + 20 format (/,' Steepest Descent Gradient Optimization :') + write (jout,30) + 30 format (/,' SD Iter F Value G RMS F Move',& + ' X Move Angle FG Call Comment',/) + else + write (jout,40) + 40 format (/,' Limited Memory BFGS Quasi-Newton',& + ' Optimization :') + write (jout,50) + 50 format (/,' QN Iter F Value G RMS F Move',& + ' X Move Angle FG Call Comment',/) + end if + flush (jout) + end if +!c +!c perform dynamic allocation of some local arrays +!c + allocate (x_old(nvar)) + allocate (g(nvar)) + allocate (g_old(nvar)) + allocate (p(nvar)) + allocate (q(nvar)) + allocate (r(nvar)) + allocate (h0(nvar)) + if (msav .ne. 0) then + allocate (rho(msav)) + allocate (alpha(msav)) + allocate (s(nvar,msav)) + allocate (y(nvar,msav)) + end if +!c +!c evaluate the function and get the initial gradient +!c + niter = nextiter - 1 + maxiter = niter + maxiter + ncalls = ncalls + 1 +! itime_mat=niter + f = fgvalue (x0,g) + f_old = f + m = 0 + gamma = 1.0d0 + g_norm = 0.0d0 + g_rms = 0.0d0 + do i = 1, nvar + g_norm = g_norm + g(i)*g(i) + g_rms = g_rms + (g(i)*scale(i))**2 + end do + g_norm = sqrt(g_norm) + g_rms = sqrt(g_rms) / rms + f_move = 0.5d0 * stpmax * g_norm +!c +!c print initial information prior to first iteration +!c + write(jout,*) "before first" + if (jprint .gt. 0) then + if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then + write (jout,60) niter,f,g_rms,ncalls + 60 format (i6,f14.4,f11.4,2x,i7) + else + write (jout,70) niter,f,g_rms,ncalls + 70 format (i6,d14.4,d11.4,2x,i7) + end if + flush (jout) + end if +!c +!c write initial intermediate prior to first iteration +!c + if (iwrite .gt. 0) call optsave (niter,f,x0) +!c +!c tests of the various termination criteria +!c + if (niter .ge. maxiter) then + statusbf = 'IterLimit' + done = .true. + end if + if (f .le. fctmin) then + statusbf = 'SmallFct ' + done = .true. + end if + if (g_rms .le. grdmin) then + statusbf = 'SmallGrad' + done = .true. + end if +!c +!c start of a new limited memory BFGS iteration +!c + do while (.not. done) + niter = niter + 1 +!c write (jout,*) "LBFGS niter",niter + muse = min(niter-1,msav) + m = m + 1 + if (m .gt. msav) m = 1 +!c +!c estimate Hessian diagonal and compute the Hg product +!c + do i = 1, nvar + h0(i) = gamma + q(i) = g(i) + end do + k = m + do j = 1, muse + k = k - 1 + if (k .eq. 0) k = msav + alpha(k) = 0.0d0 + do i = 1, nvar + alpha(k) = alpha(k) + s(i,k)*q(i) + end do + alpha(k) = alpha(k) * rho(k) + do i = 1, nvar + q(i) = q(i) - alpha(k)*y(i,k) + end do + end do + do i = 1, nvar + r(i) = h0(i) * q(i) + end do + do j = 1, muse + beta = 0.0d0 + do i = 1, nvar + beta = beta + y(i,k)*r(i) + end do + beta = beta * rho(k) + do i = 1, nvar + r(i) = r(i) + s(i,k)*(alpha(k)-beta) + end do + k = k + 1 + if (k .gt. msav) k = 1 + end do +!c +!c set search direction and store current point and gradient +!c + do i = 1, nvar + p(i) = -r(i) + x_old(i) = x0(i) + g_old(i) = g(i) + end do +!c +!c perform line search along the new conjugate direction +!c + statusbf = blank +!c write (jout,*) "Before search" + call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,statusbf) +!c write (jout,*) "After search" +!c +!c update variables based on results of this iteration +!c + if (msav .ne. 0) then + ys = 0.0d0 + yy = 0.0d0 + do i = 1, nvar + s(i,m) = x0(i) - x_old(i) + y(i,m) = g(i) - g_old(i) + ys = ys + y(i,m)*s(i,m) + yy = yy + y(i,m)*y(i,m) + end do + gamma = abs(ys/yy) + rho(m) = 1.0d0 / ys + end if +!c +!c get the sizes of the moves made during this iteration +!c + f_move = f_old - f +! if (f_move.eq.0) f_move=1.0d0 + f_old = f + x_move = 0.0d0 + do i = 1, nvar + x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2 + end do + x_move = sqrt(x_move) / rms + if (coordtype .eq. 'INTERNAL') then + x_move = radian * x_move + end if +!c +!c compute the rms gradient per optimization parameter +!c + g_rms = 0.0d0 + do i = 1, nvar + g_rms = g_rms + (g(i)*scale(i))**2 + end do + g_rms = sqrt(g_rms) / rms +!c +!c test for error due to line search problems +!c + if (statusbf.eq.'BadIntpln' .or. statusbf.eq.'IntplnErr') then + nerr = nerr + 1 + if (nerr .ge. maxerr) done = .true. + else + nerr = 0 + end if +!c +!c test for too many total iterations +!c + if (niter .ge. maxiter) then + statusbf = 'IterLimit' + done = .true. + end if +!c +!c test the normal termination criteria +!c + if (f .le. fctmin) then + statusbf = 'SmallFct ' + done = .true. + end if + if (g_rms .le. grdmin) then + statusbf = 'SmallGrad' + done = .true. + end if +!c +!c print intermediate results for the current iteration +!c +! write(iout,*) "niter1",niter + itime_mat=niter + if (jprint .gt. 0) then + if (done .or. mod(niter,jprint).eq.0) then + if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.& + g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.& + f_move.gt.-1.0d5) then + write (jout,80) niter,f,g_rms,f_move,x_move,& + angle,ncalls,statusbf + 80 format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9) + else + write (jout,90) niter,f,g_rms,f_move,x_move,& + angle,ncalls,statusbf + 90 format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9) + end if + end if + flush (jout) + end if +!c +!c write intermediate results for the current iteration +!c + if (iwrite .gt. 0) then + if (done .or. mod(niter,iwrite).eq.0) then + call optsave (niter,f,x0) + end if + end if + end do +!c +!c perform deallocation of some local arrays +!c + deallocate (x_old) + deallocate (g) + deallocate (g_old) + deallocate (p) + deallocate (q) + deallocate (r) + deallocate (h0) + if (msav .ne. 0) then + deallocate (rho) + deallocate (alpha) + deallocate (s) + deallocate (y) + end if +!c +!c set final value of the objective function +!c + minimum = f + if (jprint .gt. 0) then + if (statusbf.eq.'SmallGrad' .or. statusbf.eq.'SmallFct ') then + write (jout,100) statusbf + 100 format (/,' LBFGS -- Normal Termination due to ',a9) + else + write (jout,110) statusbf + 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9) + end if + flush (jout) + end if + return + end subroutine +!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 +! if (jjj.gt.0) then +! write (iout,*) "in func x" +! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n) +! 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) +! +! Add the components corresponding to local energy terms. +! +! 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 +!d write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg) +! g(i)=g(i)+gloc(i,icg) +! enddo +! return +! end +!endif !----------------------------------------------------------------------------- end module minimm + diff --git a/source/unres/unres.F90 b/source/unres/unres.F90 index e6c19b4..804f228 100644 --- a/source/unres/unres.F90 +++ b/source/unres/unres.F90 @@ -126,7 +126,7 @@ call exec_softreg else if (modecalc.eq.12) then call exec_MD - call exec_checkgrad +! call exec_checkgrad else if (modecalc.eq.14) then call exec_MREMD else @@ -249,7 +249,10 @@ real(kind=8) :: energy_long(0:n_ene),energy_short(0:n_ene) real(kind=8) :: varia(6*nres) !(maxvar) (maxvar=6*maxres) real(kind=8) :: time00, evals, etota, etot, time_ene, time1 - integer :: nharp,nft_sc,iretcode,nfun + integer :: nharp,nft_sc,iretcode +#ifndef LBFGS + integer :: nfun +#endif integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3) logical :: fail real(kind=8) :: rms,frac,frac_nn,co @@ -343,7 +346,13 @@ ! call check_eint ! call exec_checkgrad !el endif +! print *,'SUMSL return code is',iretcode,' eval ',nfun +#ifdef LBFGS + print *,'LBFGS return code is',status,' eval ',nfun +#else print *,'SUMSL return code is',iretcode,' eval ',nfun +#endif + #ifdef MPI evals=nfun/(MPI_WTIME()-time1) #endif @@ -382,6 +391,7 @@ use geometry_data !include 'COMMON.GEO''COMMON.CHAIN' ! use REMD !include 'COMMON.REMD' ! use MD !include 'COMMON.MD' + use regularize_ use compare implicit none @@ -409,7 +419,11 @@ if (outpdb) call pdbout(etot,titel(:32),ipdb) if (outmol2) call mol2out(etot,titel(:32)) if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) +#ifdef LBFGS + write (iout,'(a,a9)') 'LBFGS return code:',status +#else write (iout,'(a,i3)') 'SUMSL return code:',iretcode +#endif return end subroutine exec_regularize !----------------------------------------------------------------------------- @@ -494,8 +508,11 @@ real(kind=8) :: energy_(0:n_ene) logical :: eof real(kind=8) :: etot,ene0 - integer :: mm,imm,nft,n,iconf,nmin,i,iretcode,nfun,it,& + integer :: mm,imm,nft,n,iconf,nmin,i,iretcode,it,& nf_mcmf,j +#ifndef LBFGS + integer :: nfun +#endif real(kind=8) :: rms,frac,frac_nn,co,time,ene eof=.false.