adam's update
[unres.git] / source / unres / src-HCD-5D / checkder_p.F
index d1f8473..e7f0c1c 100644 (file)
@@ -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',