X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fcheckder_p.F;h=48eeddaafd2e6a532da2e8d8f73ecf8d5b98308c;hb=a30bd29e64da2aa47b84963fdd0bf4192ead2738;hp=03df2873a60242ba63ae1f35ec0238674a1940de;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;p=unres.git diff --git a/source/unres/src-HCD-5D/checkder_p.F b/source/unres/src-HCD-5D/checkder_p.F index 03df287..48eedda 100644 --- a/source/unres/src-HCD-5D/checkder_p.F +++ b/source/unres/src-HCD-5D/checkder_p.F @@ -1,181 +1,7 @@ - subroutine check_cartgrad -C Check the gradient of Cartesian coordinates in internal coordinates. - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CONTROL' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.DERIV' - dimension temp(6,maxres),xx(3),gg(3) - indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1 -* -* 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 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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. - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' @@ -183,12 +9,20 @@ C Check the gradient of the energy in Cartesian coordinates. include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTACTS' + integer i,j,k + integer icall common /srutu/ icall - dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar) - dimension grad_s(6,maxres) + double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar), + & g(maxvar),grad_s(6,maxres) 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) + double precision fdum external fdum icg=1 nf=0 @@ -202,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) @@ -253,7 +91,7 @@ C Check the gradient of the energy in Cartesian coordinates. c---------------------------------------------------------------------------- subroutine check_ecartint C Check the gradient of the energy in Cartesian coordinates. - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' @@ -264,23 +102,27 @@ C Check the gradient of the energy in Cartesian coordinates. include 'COMMON.MD' include 'COMMON.LOCAL' include 'COMMON.SPLITELE' + integer icall common /srutu/ icall - dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar), - & g(maxvar) - dimension dcnorm_safe(3),dxnorm_safe(3) - dimension grad_s(6,0:maxres),grad_s1(6,0:maxres) + 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 energia(0:n_ene),energia1(0:n_ene) integer uiparm(1) double precision urparm(1) external fdum + integer i,j,k,nf + double precision etot,etot1,etot2,etot11,etot12,etot21,etot22 + double precision dist,alpha,beta c r_cut=2.0d0 c rlambd=0.3d0 icg=1 nf=0 nfl=0 - print *,"ATU 3" +c print *,"ATU 3" call int_from_cart1(.false.) call intout c call intcartderiv @@ -325,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 @@ -349,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 @@ -493,7 +335,7 @@ c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 end c------------------------------------------------------------------------- subroutine int_from_cart1(lprn) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -509,6 +351,10 @@ c------------------------------------------------------------------------- include 'COMMON.SETUP' include 'COMMON.TIME1' logical lprn + integer i,j + double precision dnorm1,dnorm2,be + double precision time00 + double precision dist,alpha,beta if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates' #ifdef TIMING time01=MPI_Wtime() @@ -635,7 +481,7 @@ cd call flush(iout) c---------------------------------------------------------------------------- subroutine check_eint C Check the gradient of energy in internal coordinates. - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' @@ -643,14 +489,20 @@ C Check the gradient of energy in internal coordinates. include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.GEO' + integer icall common /srutu/ icall - dimension x(maxvar),gana(maxvar),gg(maxvar) + double precision x(maxvar),gana(maxvar),gg(maxvar) integer uiparm(1) double precision urparm(1) double precision energia(0:n_ene),energia1(0:n_ene), & energia2(0:n_ene) 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 c aincr=1.0D-7 print '("Calling CHECK_INT",1pd12.3)',aincr @@ -678,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 @@ -694,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',