X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fsc_move.F;h=75b7211819a628871ca0d9f3303897326a49a894;hb=c711143ad3fffb04d27b55aa823f399b8343c4c5;hp=f3535897162b2932de54df8756451544bf1bfcea;hpb=76ef494efde78d2d85d0e72d936c13166961256c;p=unres.git diff --git a/source/unres/src-HCD-5D/sc_move.F b/source/unres/src-HCD-5D/sc_move.F index f353589..75b7211 100644 --- a/source/unres/src-HCD-5D/sc_move.F +++ b/source/unres/src-HCD-5D/sc_move.F @@ -45,7 +45,7 @@ c Local variables double precision orig_w(n_ene) double precision wtime - + sideonly=.true. c Set non side-chain weights to zero (minimization is faster) c NOTE: e(2) does not actually depend on the side-chain, only CA orig_w(2)=wscp @@ -152,7 +152,8 @@ c Put the original weights back to calculate the full energy wtor=orig_w(13) wtor_d=orig_w(14) wvdwpp=orig_w(15) - + sideonly=.false. + mask_side=1 crc n_fun=n_fun+1 ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime return @@ -230,7 +231,7 @@ crc cur_e=orig_e nres_moved=0 do i=2,nres-1 c Don't do glycine (itype(j)==10) - if (itype(i).ne.10) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then sc_dist=dist(nres+i,nres+res_pick) else sc_dist=sc_dist_cutoff @@ -243,10 +244,11 @@ c Don't do glycine (itype(j)==10) endif enddo - call chainbuild + call chainbuild_extconf call egb1(evdw) call esc(escloc) e_sc=wsc*evdw+wscloc*escloc +c write (iout,*) "sc_move: e_sc",e_sc cd call etotal(energy) cd print *,'new ',(energy(k),k=0,n_ene) orig_e=e_sc @@ -271,7 +273,8 @@ crc orig_omeg(i)=omeg(i) crc enddo call minimize_sc1(e_sc,var,iretcode,loc_nfun) - +c write (iout,*) "n_try",n_try +c write (iout,*) "sc_move after minimze_sc1 e_sc",e_sc cv write(*,'(2i3,2f12.5,2i3)') cv & res_pick,nres_moved,orig_e,e_sc-cur_e, cv & iretcode,loc_nfun @@ -334,111 +337,74 @@ c Reset the minimization mask_r to false return end - -c------------------------------------------------------------- - - subroutine sc_minimize(etot,iretcode,nfun) -c Minimizes side-chains only, leaving backbone frozen -crc implicit none - -c Includes - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.FFIELD' - -c Output arguments - double precision etot - integer iretcode,nfun - -c Local variables - integer i - double precision orig_w(n_ene),energy(0:n_ene) - double precision var(maxvar) - - -c Set non side-chain weights to zero (minimization is faster) -c NOTE: e(2) does not actually depend on the side-chain, only CA - orig_w(2)=wscp - orig_w(3)=welec - orig_w(4)=wcorr - orig_w(5)=wcorr5 - orig_w(6)=wcorr6 - orig_w(7)=wel_loc - orig_w(8)=wturn3 - orig_w(9)=wturn4 - orig_w(10)=wturn6 - orig_w(11)=wang - orig_w(13)=wtor - orig_w(14)=wtor_d - - wscp=0.D0 - welec=0.D0 - wcorr=0.D0 - wcorr5=0.D0 - wcorr6=0.D0 - wel_loc=0.D0 - wturn3=0.D0 - wturn4=0.D0 - wturn6=0.D0 - wang=0.D0 - wtor=0.D0 - wtor_d=0.D0 - -c Prepare to freeze backbone - do i=1,nres - mask_phi(i)=0 - mask_theta(i)=0 - mask_side(i)=1 - enddo - -c Minimize the side-chains - mask_r=.true. - call geom_to_var(nvar,var) - call minimize(etot,var,iretcode,nfun) - call var_to_geom(nvar,var) - mask_r=.false. - -c Put the original weights back and calculate the full energy - wscp=orig_w(2) - welec=orig_w(3) - wcorr=orig_w(4) - wcorr5=orig_w(5) - wcorr6=orig_w(6) - wel_loc=orig_w(7) - wturn3=orig_w(8) - wturn4=orig_w(9) - wturn6=orig_w(10) - wang=orig_w(11) - wtor=orig_w(13) - wtor_d=orig_w(14) - - call chainbuild - call etotal(energy) - etot=energy(0) - - return - end - c------------------------------------------------------------- subroutine minimize_sc1(etot,x,iretcode,nfun) +#ifdef LBFGS_SC + use minima + use inform + use output + use iounit + use scales +#endif implicit real*8 (a-h,o-z) include 'DIMENSIONS' - parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) +#ifndef LBFGS_SC +c parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) + parameter(max_sc_move=10) + parameter (liv=60,lv=(77+2*max_sc_move*(2*max_sc_move+17)/2)) +#endif include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.GEO' include 'COMMON.MINIM' common /srutu/ icall - dimension iv(liv) - double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar) + double precision x(maxvar),d(maxvar),xx(maxvar) double precision energia(0:n_ene) +#ifdef LBFGS_SC + integer nvar_restr + common /zmienne/ nvar_restr + double precision grdmin + double precision funcgrad_restr1 + external funcgrad_restr1 + external optsave +#else external func,gradient,fdum external func_restr1,grad_restr1 logical not_done,change,reduce + dimension iv(liv) + double precision v(1:lv) common /przechowalnia/ v - +#endif +#ifdef LBFGS_SC + maxiter=7 + coordtype='RIGIDBODY' + grdmin=tolf + jout=iout +c jprint=print_min_stat + jprint=0 + 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 + call x2xx(x,xx,nvar_restr) + set_scale = .true. +c nvar = 0 + do i = 1, nvar_restr +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" + call lbfgs (nvar_restr,xx,etot,grdmin,funcgrad_restr1,optsave) + deallocate(scale) +c write (iout,*) "After lbfgs" + call xx2x(x,xx) +#else call deflt(2,iv,liv,lv,v) * 12 means fresh start, dont call deflt iv(1)=12 @@ -451,8 +417,8 @@ c------------------------------------------------------------- * controls output iv(19)=2 * selects output unit -c iv(21)=iout iv(21)=0 +c iv(21)=0 * 1 means to print out result iv(22)=0 * 1 means to print out summary stats @@ -491,14 +457,158 @@ c v(25)=4.0D0 & 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) +c call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum) ENDIF etot=v(10) iretcode=iv(1) nfun=iv(6) - +#endif return end +#ifdef LBFGS_SC + double precision function funcgrad_restr1(x,g) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.FFIELD' + include 'COMMON.INTERACT' + include 'COMMON.TIME1' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + integer nvar_restr + common /zmienne/ nvar_restr + double precision energia(0:n_ene),evdw,escloc + double precision ufparm,e1,e2 + dimension x(maxvar),g(maxvar),gg(maxvar) +#ifdef OSF +c Intercept NaNs in the coordinates, before calling etotal + x_sum=0.D0 + do i=1,nvar_restr + x_sum=x_sum+x(i) + enddo + FOUND_NAN=.false. + if (x_sum.ne.x_sum) then + write(iout,*)" *** func_restr1 : Found NaN in coordinates" + f=1.0D+73 + FOUND_NAN=.true. + return + endif +#else + FOUND_NAN=.false. + do i=1,nvar_restr + if (isnan(x(i))) then + FOUND_NAN=.true. + f=1.0D+73 + funcgrad_restr1=f + write (iout,*) "NaN in coordinates" + return + endif + enddo +#endif + +c write (iout,*) "nvar_restr",nvar_restr +c write (iout,*) "x",(x(i),i=1,nvar_restr) + call var_to_geom_restr(nvar_restr,x) + call zerograd + call chainbuild_extconf +cd write (iout,*) 'ETOTAL called from FUNC' + call egb1(evdw) + call esc(escloc) + f=wsc*evdw+wscloc*escloc +c write (iout,*) "evdw",evdw," escloc",escloc + if (isnan(f)) then + f=1.0d20 + funcgrad_restr1=f + return + endif + funcgrad_restr1=f +c write (iout,*) "f",f +cd call etotal(energia(0)) +cd f=wsc*energia(1)+wscloc*energia(12) +cd print *,f,evdw,escloc,energia(0) +C +C Sum up the components of the Cartesian gradient. +C + do i=1,nct + do j=1,3 + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i) + enddo + enddo + +C +C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +C + call cart2intgrad(nvar,gg) +C +C Convert the Cartesian gradient into internal-coordinate gradient. +C + + ig=0 + do i=2,nres-1 + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + g(ig)=gg(ialph(i,1)) +c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1) +c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)) + ENDIF + endif + enddo + do i=2,nres-1 + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + g(ig)=gg(ialph(i,1)+nside) +c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside +c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside) + ENDIF + endif + enddo + +C +C Add the components corresponding to local energy terms. +C + + ig=0 + igall=0 + do i=4,nres + igall=igall+1 + if (mask_phi(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do i=3,nres + igall=igall+1 + if (mask_theta(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + 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 + g(ig)=g(ig)+gloc(igall,icg) +c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall +c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig) + endif + endif + enddo + enddo + +cd do i=1,ig +cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +cd enddo + return + end +#else ************************************************************************ subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm) implicit real*8 (a-h,o-z) @@ -509,9 +619,7 @@ c v(25)=4.0D0 include 'COMMON.FFIELD' include 'COMMON.INTERACT' include 'COMMON.TIME1' - common /chuju/ jjj double precision energia(0:n_ene),evdw,escloc - integer jjj double precision ufparm,e1,e2 external ufparm integer uiparm(1) @@ -537,11 +645,12 @@ c Intercept NaNs in the coordinates, before calling etotal call var_to_geom_restr(n,x) call zerograd - call chainbuild + call chainbuild_extconf cd write (iout,*) 'ETOTAL called from FUNC' call egb1(evdw) call esc(escloc) f=wsc*evdw+wscloc*escloc +c write (iout,*) "f",f cd call etotal(energia(0)) cd f=wsc*energia(1)+wscloc*energia(12) cd print *,f,evdw,escloc,energia(0) @@ -550,7 +659,7 @@ C Sum up the components of the Cartesian gradient. C do i=1,nct do j=1,3 - gradx(j,i,icg)=wsc*gvdwx(j,i) + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i) enddo enddo @@ -569,7 +678,7 @@ c------------------------------------------------------- external ufparm integer uiparm(1) double precision urparm(1) - dimension x(maxvar),g(maxvar) + dimension x(maxvar),g(maxvar),gg(maxvar) icg=mod(nf,2)+1 if (nf-nfl+1) 20,30,40 @@ -578,76 +687,51 @@ c write (iout,*) 'grad 20' if (nf.eq.0) return goto 40 30 call var_to_geom_restr(n,x) - call chainbuild + call chainbuild_extconf C C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. C - 40 call cartder + 40 call cart2intgrad(nvar,gg) C C Convert the Cartesian gradient into internal-coordinate gradient. C ig=0 - ind=nres-2 + ind=nres-2 do i=2,nres-2 - IF (mask_phi(i+2).eq.1) THEN - gphii=0.0D0 - do j=i+1,nres-1 - ind=ind+1 - do k=1,3 - gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) - gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) - enddo - enddo + IF (mask_phi(i+2).eq.1) THEN ig=ig+1 - g(ig)=gphii - ELSE - ind=ind+nres-1-i + g(ig)=gg(i-1) ENDIF enddo - ind=0 do i=1,nres-2 IF (mask_theta(i+2).eq.1) THEN ig=ig+1 - gthetai=0.0D0 - do j=i+1,nres-1 - ind=ind+1 - do k=1,3 - gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) - gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) - enddo - enddo - g(ig)=gthetai - ELSE - ind=ind+nres-1-i + g(ig)=gg(nphi+i) ENDIF enddo do i=2,nres-1 - if (itype(i).ne.10) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then IF (mask_side(i).eq.1) THEN ig=ig+1 - galphai=0.0D0 - do k=1,3 - galphai=galphai+dxds(k,i)*gradx(k,i,icg) - enddo - g(ig)=galphai + g(ig)=gg(ialph(i,1)) +c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1) +c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)) ENDIF endif enddo do i=2,nres-1 - if (itype(i).ne.10) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then IF (mask_side(i).eq.1) THEN ig=ig+1 - gomegai=0.0D0 - do k=1,3 - gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) - enddo - g(ig)=gomegai + g(ig)=gg(ialph(i,1)+nside) +c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside +c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside) ENDIF endif enddo @@ -676,11 +760,13 @@ C do ij=1,2 do i=2,nres-1 - if (itype(i).ne.10) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then igall=igall+1 if (mask_side(i).eq.1) then ig=ig+1 g(ig)=g(ig)+gloc(igall,icg) +c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall +c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig) endif endif enddo @@ -691,6 +777,7 @@ cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) cd enddo return end +#endif C----------------------------------------------------------------------------- subroutine egb1(evdw) C @@ -716,11 +803,12 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do i=nnt,nct itypi=iabs(itype(i)) - if (itypi.eq.ntyp1) cycle + if (itypi.eq.ntyp1 .or. mask_side(i).eq.0) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -761,8 +849,9 @@ C lipbufthick is thickenes of lipid buffore C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) + do j=i+1,nct IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN ind=ind+1 itypj=iabs(itype(j)) @@ -922,7 +1011,7 @@ C Calculate angular part of the gradient. call sc_grad ENDIF enddo ! j - enddo ! iint +c enddo ! iint enddo ! i end C-----------------------------------------------------------------------------