X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc-HCD-5D%2Fcheckder_p.F;h=48eeddaafd2e6a532da2e8d8f73ecf8d5b98308c;hb=8c073a2ffd0a50442e7fdf20780a0dd94ed74799;hp=d1f84730cd101ec20012700faf5919a38eb8e309;hpb=48ae9e01d2dd6571fa2cca6c704dc04f86e5fd7b;p=unres.git diff --git a/source/unres/src-HCD-5D/checkder_p.F b/source/unres/src-HCD-5D/checkder_p.F index d1f8473..48eedda 100644 --- a/source/unres/src-HCD-5D/checkder_p.F +++ b/source/unres/src-HCD-5D/checkder_p.F @@ -1,182 +1,3 @@ - subroutine check_cartgrad -C Check the gradient of Cartesian coordinates in internal coordinates. - implicit none - include 'DIMENSIONS' - include 'COMMON.CONTROL' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.DERIV' - double precision temp(6,maxres),xx(3),gg(3),thet,theti,phii,alphi, - & omegi,aincr2 - integer indmat - integer i,ii,j,k - indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1 - integer nf -* -* Check the gradient of the virtual-bond and SC vectors in the internal -* coordinates. -* - print '("Calling CHECK_ECART",1pd12.3)',aincr - write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr - aincr2=0.5d0*aincr - call chainbuild_extconf - call cartder - write (iout,'(a)') '**************** dx/dalpha' - write (iout,'(a)') - do i=2,nres-1 - alphi=alph(i) - alph(i)=alph(i)+aincr - do k=1,3 - temp(k,i)=dc(k,nres+i) - enddo - call chainbuild_extconf - do k=1,3 - gg(k)=(dc(k,nres+i)-temp(k,i))/aincr - xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr)) - enddo - write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') - & i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3) - write (iout,'(a)') - alph(i)=alphi - call chainbuild_extconf - enddo - write (iout,'(a)') - write (iout,'(a)') '**************** dx/domega' - write (iout,'(a)') - do i=2,nres-1 - omegi=omeg(i) - omeg(i)=omeg(i)+aincr - do k=1,3 - temp(k,i)=dc(k,nres+i) - enddo - call chainbuild_extconf - do k=1,3 - gg(k)=(dc(k,nres+i)-temp(k,i))/aincr - xx(k)=dabs((gg(k)-dxds(k+3,i))/ - & (aincr*dabs(dxds(k+3,i))+aincr)) - enddo - write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') - & i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3) - write (iout,'(a)') - omeg(i)=omegi - call chainbuild_extconf - enddo - write (iout,'(a)') - write (iout,'(a)') '**************** dx/dtheta' - write (iout,'(a)') - do i=3,nres - theti=theta(i) - theta(i)=theta(i)+aincr - do j=i-1,nres-1 - do k=1,3 - temp(k,j)=dc(k,nres+j) - enddo - enddo - call chainbuild_extconf - do j=i-1,nres-1 - ii = indmat(i-2,j) -c print *,'i=',i-2,' j=',j-1,' ii=',ii - do k=1,3 - gg(k)=(dc(k,nres+j)-temp(k,j))/aincr - xx(k)=dabs((gg(k)-dxdv(k,ii))/ - & (aincr*dabs(dxdv(k,ii))+aincr)) - enddo - write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') - & i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3) - write(iout,'(a)') - enddo - write (iout,'(a)') - theta(i)=theti - call chainbuild_extconf - enddo - write (iout,'(a)') '***************** dx/dphi' - write (iout,'(a)') - do i=4,nres - phi(i)=phi(i)+aincr - do j=i-1,nres-1 - do k=1,3 - temp(k,j)=dc(k,nres+j) - enddo - enddo - call chainbuild_extconf - do j=i-1,nres-1 - ii = indmat(i-2,j) -c print *,'ii=',ii - do k=1,3 - gg(k)=(dc(k,nres+j)-temp(k,j))/aincr - xx(k)=dabs((gg(k)-dxdv(k+3,ii))/ - & (aincr*dabs(dxdv(k+3,ii))+aincr)) - enddo - write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') - & i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3) - write(iout,'(a)') - enddo - phi(i)=phi(i)-aincr - call chainbuild_extconf - enddo - write (iout,'(a)') '****************** ddc/dtheta' - do i=1,nres-2 - thet=theta(i+2) - theta(i+2)=thet+aincr - do j=i,nres - do k=1,3 - temp(k,j)=dc(k,j) - enddo - enddo - call chainbuild_extconf - do j=i+1,nres-1 - ii = indmat(i,j) -c print *,'ii=',ii - do k=1,3 - gg(k)=(dc(k,j)-temp(k,j))/aincr - xx(k)=dabs((gg(k)-dcdv(k,ii))/ - & (aincr*dabs(dcdv(k,ii))+aincr)) - enddo - write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') - & i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3) - write (iout,'(a)') - enddo - do j=1,nres - do k=1,3 - dc(k,j)=temp(k,j) - enddo - enddo - theta(i+2)=thet - enddo - write (iout,'(a)') '******************* ddc/dphi' - do i=1,nres-3 - phii=phi(i+3) - phi(i+3)=phii+aincr - do j=1,nres - do k=1,3 - temp(k,j)=dc(k,j) - enddo - enddo - call chainbuild_extconf - do j=i+2,nres-1 - ii = indmat(i+1,j) -c print *,'ii=',ii - do k=1,3 - gg(k)=(dc(k,j)-temp(k,j))/aincr - xx(k)=dabs((gg(k)-dcdv(k+3,ii))/ - & (aincr*dabs(dcdv(k+3,ii))+aincr)) - enddo - write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') - & i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3) - write (iout,'(a)') - enddo - do j=1,nres - do k=1,3 - dc(k,j)=temp(k,j) - enddo - enddo - phi(i+3)=phii - enddo - return - end C---------------------------------------------------------------------------- subroutine check_ecart C Check the gradient of the energy in Cartesian coordinates. @@ -196,6 +17,8 @@ C Check the gradient of the energy in Cartesian coordinates. double precision energia(0:n_ene),energia1(0:n_ene) double precision aincr2,etot,etot1,etot2 double precision dist,alpha,beta + double precision funcgrad,ff + external funcgrad integer nf integer uiparm(1) double precision urparm(1) @@ -213,7 +36,11 @@ C Check the gradient of the energy in Cartesian coordinates. call etotal(energia(0)) etot=energia(0) call enerprint(energia(0)) +#ifdef LBFGS + ff=funcgrad(x,g) +#else 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) @@ -261,259 +88,6 @@ C Check the gradient of the energy in Cartesian coordinates. enddo return end -#ifdef FIVEDIAG -!----------------------------------------------------------------------------- - subroutine check_ecartint -! Check the gradient of the energy in Cartesian coordinates. - implicit none - include 'DIMENSIONS' - include 'COMMON.CONTROL' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.DERIV' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CONTACTS' - include 'COMMON.MD' - include 'COMMON.LOCAL' - include 'COMMON.SPLITELE' - integer icall - common /srutu/ icall - double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3), - & x(maxvar),g(maxvar) - double precision dcnorm_safe(3),dxnorm_safe(3) - double precision grad_s(6,0:maxres),grad_s1(6,0:maxres) - double precision phi_temp(maxres),theta_temp(maxres), - & alph_temp(maxres),omeg_temp(maxres) - double precision ddc1(3),ddcn(3),dcnorm_safe1(3),dcnorm_safe2(3) - double precision energia(0:n_ene),energia1(0:n_ene) - integer uiparm(1) - double precision urparm(1) - double precision fdum - external fdum - integer i,j,k,nf - double precision etot,etot1,etot2,etot11,etot12,etot21,etot22 - double precision dist,alpha,beta - icg=1 - nf=0 - nfl=0 - call intout -! call intcartderiv -! call checkintcartgrad - call zerograd - aincr=1.0D-5 - write(iout,*) 'Calling CHECK_ECARTINT.' - nf=0 - icall=0 - write (iout,*) "Before geom_to_var" - call geom_to_var(nvar,x) - write (iout,*) "after geom_to_var" - write (iout,*) "split_ene ",split_ene - call flush(iout) - if (.not.split_ene) then - write(iout,*) 'Calling CHECK_ECARTINT if' - call etotal(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - etot=energia(0) - write (iout,*) "etot",etot - call enerprint(energia(0)) - call flush(iout) -!el call enerprint(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) - call cartgrad -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - write (iout,*) "exit cartgrad" - call flush(iout) - icall =1 - do i=1,nres - write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) - enddo - do j=1,3 - grad_s(j,0)=gcart(j,0) - enddo -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - do i=1,nres - do j=1,3 - grad_s(j,i)=gcart(j,i) - grad_s(j+3,i)=gxcart(j,i) - enddo - enddo - else - write(iout,*) 'Calling CHECK_ECARTIN else.' -!- split gradient check - call zerograd - call etotal_long(energia) - call enerprint(energia(0)) -!el call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) - call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) - icall =1 - write (iout,*) "longrange grad" - do i=1,nres - write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), - & (gxcart(j,i),j=1,3) - enddo - do j=1,3 - grad_s(j,0)=gcart(j,0) - enddo - do i=1,nres - do j=1,3 - grad_s(j,i)=gcart(j,i) - grad_s(j+3,i)=gxcart(j,i) - enddo - enddo - call zerograd - call etotal_short(energia) - call enerprint(energia(0)) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) - call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) - icall =1 - write (iout,*) "shortrange grad" - do i=1,nres - write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), - & (gxcart(j,i),j=1,3) - enddo - do j=1,3 - grad_s1(j,0)=gcart(j,0) - enddo - do i=1,nres - do j=1,3 - grad_s1(j,i)=gcart(j,i) - grad_s1(j+3,i)=gxcart(j,i) - enddo - enddo - endif - write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' -! do i=1,nres -c do i=nnt,nct - do i=1,nres - do j=1,3 - if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1) - if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres) - ddc(j)=c(j,i) - ddx(j)=c(j,i+nres) - dcnorm_safe1(j)=dc_norm(j,i-1) - dcnorm_safe2(j)=dc_norm(j,i) - dxnorm_safe(j)=dc_norm(j,i+nres) - enddo - do j=1,3 - c(j,i)=ddc(j)+aincr - if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr - if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr - if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) - dc(j,i)=c(j,i+1)-c(j,i) - dc(j,i+nres)=c(j,i+nres)-c(j,i) - call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1) - etot1=energia1(0) -c write (iout,*) "ij",i,j," etot1",etot1 - else -!- split gradient - call etotal_long(energia1) - etot11=energia1(0) - call etotal_short(energia1) - etot12=energia1(0) - endif -!- end split gradient -! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 - c(j,i)=ddc(j)-aincr - if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr - if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr - if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) - dc(j,i)=c(j,i+1)-c(j,i) - dc(j,i+nres)=c(j,i+nres)-c(j,i) - call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1) - etot2=energia1(0) -c write (iout,*) "ij",i,j," etot2",etot2 - ggg(j)=(etot1-etot2)/(2*aincr) - else -!- split gradient - call etotal_long(energia1) - etot21=energia1(0) - ggg(j)=(etot11-etot21)/(2*aincr) - call etotal_short(energia1) - etot22=energia1(0) - ggg1(j)=(etot12-etot22)/(2*aincr) -!- end split gradient -! write (iout,*) "etot21",etot21," etot22",etot22 - endif -! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 - c(j,i)=ddc(j) - if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j) - if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j) - if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) - dc(j,i)=c(j,i+1)-c(j,i) - dc(j,i+nres)=c(j,i+nres)-c(j,i) - dc_norm(j,i-1)=dcnorm_safe1(j) - dc_norm(j,i)=dcnorm_safe2(j) - dc_norm(j,i+nres)=dxnorm_safe(j) - enddo - do j=1,3 - c(j,i+nres)=ddx(j)+aincr - dc(j,i+nres)=c(j,i+nres)-c(j,i) - call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1) - etot1=energia1(0) - else -!- split gradient - call etotal_long(energia1) - etot11=energia1(0) - call etotal_short(energia1) - etot12=energia1(0) - endif -!- end split gradient - c(j,i+nres)=ddx(j)-aincr - dc(j,i+nres)=c(j,i+nres)-c(j,i) - call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1) - etot2=energia1(0) - ggg(j+3)=(etot1-etot2)/(2*aincr) - else -!- split gradient - call etotal_long(energia1) - etot21=energia1(0) - ggg(j+3)=(etot11-etot21)/(2*aincr) - call etotal_short(energia1) - etot22=energia1(0) - ggg1(j+3)=(etot12-etot22)/(2*aincr) -!- end split gradient - endif -! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 - c(j,i+nres)=ddx(j) - dc(j,i+nres)=c(j,i+nres)-c(j,i) - dc_norm(j,i+nres)=dxnorm_safe(j) - call int_from_cart1(.false.) - enddo - write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') - & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6) - if (split_ene) then - write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') - & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i), - & k=1,6) - write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') - & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6), - & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6) - endif - enddo - return - end -#else c---------------------------------------------------------------------------- subroutine check_ecartint C Check the gradient of the energy in Cartesian coordinates. @@ -593,15 +167,15 @@ c call flush(iout) call etotal_long(energia(0)) call enerprint(energia(0)) call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) +c write (iout,*) "enter cartgrad" +c call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) +c write (iout,*) "exit cartgrad" +c call flush(iout) icall =1 write (iout,*) "longrange grad" do i=1,nres - write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), + write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo do j=1,3 @@ -617,15 +191,15 @@ c call flush(iout) call etotal_short(energia(0)) call enerprint(energia(0)) call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) +c write (iout,*) "enter cartgrad" +c call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) +c write (iout,*) "exit cartgrad" +c call flush(iout) icall =1 write (iout,*) "shortrange grad" do i=1,nres - write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), + write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo do j=1,3 @@ -759,7 +333,6 @@ c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 enddo return end -#endif c------------------------------------------------------------------------- subroutine int_from_cart1(lprn) implicit none @@ -926,6 +499,8 @@ C Check the gradient of energy in internal coordinates. character*6 key double precision fdum external fdum + double precision funcgrad,ff + external funcgrad integer i,ii,nf double precision xi,etot,etot1,etot2 call zerograd @@ -955,7 +530,15 @@ c aincr=1.0D-7 nf=1 nfl=3 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) +c write (iout,*) "Before gradient" +c call flush(iout) +#ifdef LBFGS + ff=funcgrad(x,gana) +#else call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) +#endif +c write (iout,*) "After gradient" +c call flush(iout) cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) icall=1 do i=1,nvar @@ -971,7 +554,7 @@ cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) call etotal(energia2(0)) etot2=energia2(0) gg(i)=(etot2-etot1)/aincr - write (iout,*) i,etot1,etot2 +c write (iout,*) i,etot1,etot2 x(i)=xi enddo write (iout,'(/2a)')' Variable Numerical Analytical',