X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fcheckder_p.F;h=e7f0c1c772a11c2f7752cf6938e0ef1424e8756b;hb=c711143ad3fffb04d27b55aa823f399b8343c4c5;hp=d1f84730cd101ec20012700faf5919a38eb8e309;hpb=76ef494efde78d2d85d0e72d936c13166961256c;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..e7f0c1c 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) @@ -320,17 +147,19 @@ C Check the gradient of the energy in Cartesian coordinates. call flush(iout) !el call enerprint(energia) !elwrite(iout,*) 'Calling CHECK_ECARTINT if' - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) +c call flush(iout) +c write (iout,*) "enter cartgrad" +c call flush(iout) call cartgrad +c Transform the gradient to the CA-SC basis + call grad_transform !elwrite(iout,*) 'Calling CHECK_ECARTINT if' - write (iout,*) "exit cartgrad" - call flush(iout) +c write (iout,*) "exit cartgrad" +c call flush(iout) icall =1 - do i=1,nres - write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) - enddo +c do i=1,nres +c write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) +c enddo do j=1,3 grad_s(j,0)=gcart(j,0) enddo @@ -342,18 +171,20 @@ C Check the gradient of the energy in Cartesian coordinates. enddo enddo else - write(iout,*) 'Calling CHECK_ECARTIN else.' +c 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) +c call flush(iout) +c write (iout,*) "enter cartgrad" +c call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) +c Transform the gradient to CA-SC coordinates + call grad_transform +c write (iout,*) "exit cartgrad" +c call flush(iout) icall =1 write (iout,*) "longrange grad" do i=1,nres @@ -373,11 +204,13 @@ C Check the gradient of the energy in Cartesian coordinates. call etotal_short(energia) 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) +c Transform the gradient to CA-SC basis + call grad_transform icall =1 write (iout,*) "shortrange grad" do i=1,nres @@ -926,6 +759,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 +790,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 +814,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',